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■'  A  series  of  computer  programs  written  in  ICL  1900  series  FORTRAN  is 
presented  to  generate  data  for  the  finite  element  analysis  of  shells,  with  parti¬ 
cular  reference  to  spherical  surfaces.  The  programs  are  appropriate  to  a  version 
of  the  SEMILOOF  element  contained  in  an  RAE  Structures  Department  program  and  to 
TRIA3  and  QUAD4  elements  which  are  available  in  a  NASTRAN  package.  Both 
descriptions  and  listings  of  the  programs  are  given.  ~ 
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1  INTRODUCTION 

The  data  input  to  finite  element  packages  is  often  of  a  repetitive  nature 
and  in  fixed  format  which  makes  manual  data  preparation  both  time-consuming  and 
prone  to  error.  The  computer  programs  described  and  listed  here  form  a  less 
tedious  and  faster  alternative  to  this  manual  stage  of  data  preparation. 

The  programs  are  appropriate  either  to  a  finite  element  program  containing 
the  SEMILOOF  thin  shell  element,  held  by  Structures  Department,  RAE  (see 
the  Appendix)  or  for  a  NASTRAN  package*,  containing  the  TRIA3  and  QUAD4  facet 
shell  elements,  which  are  accessed  via  a  commercial  bureau.  The  majority  of  the 
programs  generate  data  specifically  for  meshes  on  the  surface  of  a  sphere,  and 

they  have  been  used  to  generate  all  the  data  for  the  hemispherical  cap  and 

.  2 
patch  test  problems  considered  by  Morley  and  Morris  . 

Two  types  of  program  are  described.  One  type  alters  a  free  format  data 
file  to  one  that  can  be  incorporated  directly  into  a  data  deck.  The 
other  type  is  for  a  particular  mesh,  ie  a  ’patch',  in  which  only  the  essential 
variable  parameters  are  input. 

2  PROGRAM  GENDATA 

This  program  generates  the  x,  y  and  z  coordinates  on  the  surface  of  a 
sphere  of  radius  R  when  the  corresponding  0  and  ij>  surface  coordinates  are 
given.  The  program  can  be  used  to  create  either  SEMILOOF  or  NASTRAN  QUAD4  data. 
In  the  latter  case  the  GENDATA  output  is  used  as  input  to  program  GENGRID.  For 
certain  SEMILOOF  meshes  it  is  convenient  to  interpolate  the  midside  node 
positions  from  vertex  nodes.  For  this  option,  dummy  values  of  0  and  <p  at 
the  midside  nodes  are  given  initially  and,  when  the  calculation  of  the  rest  of 
the  nodes  is  completed,  the  positions  of  these  midside  nodes  are  calculated. 

This  is  done  by  reading  in  the  number  of  the  midside  node  together  with  the  nodes 
between  which  it  lies  and  bisecting  the  geodesic  line  between  the  vertices.  This 
is  repeated  until  the  end  of  the  file  is  reached.  Note  that  constant  0  gives  a 
"line  of  longitude"  which  is  a  geodesic  line,  but  constant  i$>  gives  a  "line  of 
latitude"  which  is  not  a  geodesic  line. 

The  radius  of  the  sphere,  R  ,  is  set  to  10. 0  but  this  can  be  altered  by 
changing  a  single  line  in  the  program. 
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The  free  format  input  data  is  as  followsi- 

NNODE 

NODE  |  6j  $ 

•  •  • 


NODE.  6. 

1  11 


NODE 

nnode 

NODE 


where  NNODE  is  the  total  number  of  vertex  and  midside  nodes, 

0;  and  <f^  are  the  6  and  4>  coordinates  in  degrees  of  node  i,  1,2,..., 
i, . . . , NNODE, 

NODE  is  the  midside  node  lying  between  NOA  and  NOB. 

The  output  is  to  card  punch  CPO  and  is  of  the  form  NODE^ 
where  X,  Y,  Z  are  the  x  y  z  coordinates  of  NODE^  1 ,2, . . . ,i, . . . .NNODE. 

The  nodes  are  input  in  sequential  order  and  checked.  If  a  node  is  found  to  be 
out  of  sequence  a  message  is  output  to  line  printer  LPO. 

3  PATCH  DATA  GENERATION  PROGRAMS 

These  programs  generate  the  x,  y  and  z  coordinates,  and  boundary  condition 
data  for  a  patch  of  nine  quadrilateral  elements  or  18  triangular  elements  on  the 
surface  of  a  sphere  of  radius  R  -  10.0,  see  Fig  4. 


Programs  GENPATCH  and  GENTRIPAT 


Program  GENPATCH  calculates  the  data  for  a  3  x  3  regular  mesh  of  SEMILOOF 
or  NASTRAN  QUAD4  elements  as  shown  in  Fig  2a&b  respectively.  Program  GENTRIPAT 
calculates  the  data  for  a  similar  mesh  of  triangular  SEMILOOF  elements  as  in 
Fig  3a.  The  input  for  both  of  these  programs  is  the  same,  ie 


ITYPE  ICASE  THE INC  PHI INC  THEZERO  PHIZERO 


where  ITYPE  »  1  for  SEMILOOF  output, 


2  for  NASTRAN  QUAD4  output, 
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ICASE  is  a  variable  used  in  subsequent  programs  to  specify  the  type  of 
boundary  displacement  to  be  applied, 

»  0  for  pressure  loading, 

=  1  for  k  constant, 

■  2  for  k^,  -Kq9  constant, 

=  3  for  e  constant. 


=  3  for 
=  4  for 


constant. 


are  the  mean  of  the  values  of  6  and 


THEINC  is  the  increment  per  element  in  the  0  direction,  of  Fig  4b, 

PHIINC  is  the  increment  per  element  in  the  $  direction*  of  Fig  4b, 

THEZERO  is  the  6  coordinate  of  the  centre  of  the  patch  (usually  0.0  or 

45.0), 

PH1ZER0  is  the  <t>  coordinate  of  the  centre  of  the  patch  (usually  90.0). 

THEINC,  PHIINC,  THEZERO  and  PHIZER0  are  given  in  degrees. 

During  the  execution  of  these  programs  arrays  X,  Y,  Z  are  filled  with  the 
x,  y  and  z  coordinates  of  each  node.  Arrays  ATHE  and  APHI  are  filled  with 
the  6  and  $  coordinates  in  degrees.  Note  that  the  midside  values  of  0  and 
<j>  are  the  mean  of  the  values  of  0  and  <|>  at  the  vertex  nodes,  thus  the  mid¬ 
side  node  is  not  necessarily  on  the  geodesic  line  between  the  vertex  nodes. 

For  SEMILOOF  data,  ie  ITYPE  =  1 ,  the  output  to  CPO  is 


for  i  *  1,  2,...i,...,  NNODE .  The  output  to  CPI  is 


ICASE 


M  0„  I  SIDE 

where  N  is  repeated  for  each  vertex  node  around  the  boundary  and  M  is 
repeated  for  each  midside  node  around  the  boundary.  ISIDE  =  1,  2,  3  or  4  accord¬ 
ing  to  the  side  (see  Fig  3a) . 

For  NASTRAN  QUAD4  data,  ie  ITYPE  =  2,  the  output  to  CPO  is 

GRID  NODE.  X.  Y.  Z. 

1  1  1.1 

for  i  ■  1,  2,...,  i,...,  NNODE.  This  can  be  incorporated  directly  into  a  NASTRAN 
bulk  data  deck.  The  output  to  CPI  is 


ICASE 


ISIDE 


where  N  is  repeated  for  each  node  on  the  boundary  and  ISIDE  =  1,  2,  3,  4  accord¬ 
ing  to  the  side  and  ISIDE  is  set  equal  to  5  for  a  node  on  the  corner  of  the  patch. 
This  can  then  be  used  directly  as  input  to  program  GENNASTBC  to  obtain  the  actual 
displacements  around  the  boundary. 


3.2  Program  GBNIPAT 

This  program  generates  data  for  a  3  *  3  patch  of  SEMILOOF  or  NASTRAN  QUAD4 
elements  on  the  surface  of  a  sphere.  The  8  and  $  coordinates  of  the  four 
central  vertex  nodes  are  given,  enabling  the  user  to  generate  the  coordinates  of 
a  mesh  with  a  regular  boundary  but  an  irregular  interior,  eg  as  in  Fig  3c. 

The  first  line  of  input  is  identical  to  the  input  for  program  GENPATCH, 
as  described  in  section  3.1,  ie 

ITYPE  ICASE  THEINC  PHIINC  THEZERO  PHIZERO 
the  subsequent  lines  are,  for  SEMILOOF: 

0, 


14 


14 


0 


16 


♦ 


16 


25 


e 


27 


25 

27 


and  for  NASTRAN: 

0. 


0. 


10 


0 


11 


♦ 7 

*10 

*11 


All  0  and  4  coordinates  are  given  in  degrees. 

The  output  is  similar  to  that  of  program  GENPATCH. 


4  PROGRAM  GENNASTBC 

This  program  generates  the  prescribed  boundary  condition  data  for  a  3  *  3 
patch  of  NASTRAN  QUAD4  elements  on  the  surface  of  a  sphere  as  in  Fig  2b.  This 
can  be  done  for  each  of  five  patch  cases  as  specified  by  an  input  parameter. 

The  input  to  the  program  is  the  output  on  CPI  from  program  GENPATCH,  ie 


ICASE 

N  0N  4N  I  SIDE 


where  ICASE  » 

0  for 

constant  pressure  loading. 

■ 

1  for 

%  constant, 

■ 

2  for 

*44’  "*00  con8tant’ 

m 

3  for 

e0^  constant, 

m 

4  for 

e44  *  "eee  constant, 

N 

is  the 

node  number  on  the  boundary, 

6*,  and  4^  are  the  0  and  $  coordinates,  in  degrees  at  node  N, 

I SIDE  is  the  side  number  1,  2,  3  or  4  on  which  the  node  is  situated 

or  5  if  the  node  is  on  the  corner  of  a  patch,  (ISIDE  is  only 
read  for  ICASE  ■  0) . 
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The  output  is  to  three  card  punch  files  CPO,  CPI,  CP2  ready  for  inclusion 
in  a  NASTRAN  bulk  data  deck. 

4.1  Subroutine  PRESSURE 

This  subroutine  calculates  the  boundary  conditions  for  a  patch  in  a  state 
of  constant  pressure  acting  in  a  complete  sphere.  The  subroutine  is  called  when 
1CASE  =  0.  For  unit  pressure  loading  U,  ,  U  and  W  ,  the  displacements  in  the 

(f)  U 

<t>  and  6  directions  and  outward  normal  to  the  shell  surface  respectively,  are 


=  0 


W 


=  0  , 


kR 


2D  ( 1  +  v)  J 


(4-1) 


*  •! 

SCJ 


f 


where  k  = 


I2R 


h  is  the  shell  thickness 
D  is  the  flexural  rigidity 

Eh3 


12(1  -  v  ) 

v  is  Poisson's  ratio. 

The  patch  may  be  centred  anywhere  on  the  surface  of  the  sphere.  Conditions 
of  symmetry  are  applied  to  each  side  of  the  patch  and  are  output  on  single 
precision  MPC  (multi-point  constraint)  cards  on  card  punch  file  CPO.  (Files  CPI 
and  CP2  remain  empty.)  On  sides  1  and  3  (see  Fig  3a)  and  5,  (where  5  is  a 
corner  of  the  patch)  the  conditions  are 


-  cos  <4  cos  e  u  -  cos  it  sin  6  U  +  sin  <4  U 
x  y  T  z 


-  sin  0  w  +  cos  0  w 
x  y 


=  0 


(4-2) 


and  on  sides  2,  4  and  5  they  are 


-  sin  0  U  +  cos  0  U  “0 
x  y 


cos  <4  cos  0  w  -  cos  di  sine  u  +  sin  <4  w  ■  0 

x  T  y  z 


(4-3) 
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where  U  ,  U  ,  U  and  to  ,  tu  ,  a)  are  the  displacements  and  rotations  in  the 
x  y  z  x  y  z  r 

x,  y  and  z  directions. 

4.2  Subroutine  KTHEPHI 

This  subroutine  is  called  when  ICASE  "I.  It  outputs  the  prescribed 
inextensional  bending  case  displacements  for  the  nearly  constant  curvature  change 
state,  where 

—  constant  ^ 


(4-4) 


K  =  -  <  —  0  . 

<t»b  00 


The  prescribed  displacements  are 


(4-5) 


=  y  sin  <f>  ^tan2  +  cot2  ^  cos  20  , 


*  B(2  -  sin  ij>)  cosec  $  cos  20  , 


U„  *  -  2B  cot  <J>  sin  20  , 


W  <=  B(2  +  sin  d>)  cot  4>  cosec  cos  20 


When  the  patch  is  centred  at  0  *  45°,  <f>  =  90°  these  displacements  may  be  written 


3 

15  4B  cos  <(>  cos  0/sin  <f>  , 


3 

Uy  *  -  4B  cos  <j>  sin  0/sin  <f>  , 


2  2  2 

U  **  -  B  cos  2  0  (sin  d>  —  2  cos  4>)/sin  41 

z 


(4-6) 


a)  *  4B  cos2  0/R  sin  41  , 


*  4B  sin  2  0  cos  <f>/R  sin  4  , 


where  B  =  R2(|  -  v)/Eh^ 

and  w  and  w  are  the  rotations  about  the  0  and  $  directions. 
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The  quantities  U  ,  U  and  U  ,  ie  freedoms  1,  2  and  3  are  prescribed  on 
x  y  z 

double  precision  SPC  cards  on  CPO.  id.  and  u»  are  resolved  into  w  ,  w  and 

0  <J>  x  y 

(oz  and  are  prescribed  on  double  precision  MPC  cards  on  CPI,  ie 


-  sin  6  to  +  cos  0ui  -  o),  =  0  , 

x  y  <t> 


-  cos  0  cos  $  o)  -  cos  6  sin0  id  +  sin<(i  id  -  id.  =  0  . 

x  y  z  0 


(4-7) 


To  represent  the  constant  terms  oj  and  id  in  the  multipoint  constraints  it 

0  9 

is  necessary  to  introduce  scalar  points  which  are  prescribed  to  have  unit  dis¬ 
placement  on  single  precision  SPC  cards  on  CP2. 

4.3  Subroutine  KTHETHE 

This  subroutine,  called  when  ICASE  =  2  ,  is  for  the  nearly  constant  curva¬ 
ture  change  state 


where 


K  =  -  K 


K  -  0  . 


constant  , 


(4-8) 


The  prescribed  displacements  are 


=  ~2  sin  $  ^tan2  -  cot2  cos  20  , 


=  -  —  cot  <(>  cos  20  , 


(4-9) 


1  2 

U„  =  -r  (2  -  sin  <t)  cosec  $  sin  20  , 
o  o 

1  2 

W  =  -  —  cosec  q>  cos  20  . 

i  i 

When  the  patch  is  centred  at  0  =  0°,  $  =  90°  these  displacements  may  be  written 


3  2 

-  cos  0(cos  <f>  +  l)/3  sin  <j>  , 
3  2 

sin  0(cos  +  l)/3  sin  cp  , 


cos  20  cos  <f> 


(' '  ~t)/ 

\  sin  <py 


(4-10) 
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2  3 

u)^  *  -  cos  <j>  cos  28  (sin  <j>  +  2)/3R  sin  $  , 

.24  3 

a)  *  sin  20(2  sin  <J>  -  sin  <J>  -  4)/6R  sin  <j>  . 

U 

The  quantities  U^,  U^,  are  prescribed  on  double  precision  SPC  cards  on  CPO. 

m  and  w.  are  resolved,  as  in  equation  (4-7),  into  o>  ,  w  and  oi  and  are 
0  9  x  y  z 

prescribed  on  double  precision  MPC  cards  on  CPI.  As  in  subroutine  KTHEPHI 
scalar  points  are  again  introduced  and  their  displacements  fixed  on  single 
precision  SPC  cards  output  to  CP2. 

4.4  Subroutine  ETHEPHI 

This  subroutine,  called  for  ICASE  =  3,  prescribes  boundary  displacements 
for  the  nearly  constant  strain  patch  test 


e..  ~  constant 


e, .  *  -  e. 


(4-11) 


The  appropriate  displacements  are 


=  -  ~  £|3  +  (2  +  cos  <|>)  cos  <(>|  tan2  -| 

+  ^3  -  (2  -  cos  <j>)  cos  <j>  |  cot2  -|-J  cosec  $  cos  20  , 

13  .2  4 

=  -  -g  cosec  (J>(4  -  2  sin  $  +  sin  <f>)  cos  20  , 


(4-12) 


1  2  2 
U0  =  -  j  cot  <t>  cosec  ij>(2  +  sin  <p)  sin  20  , 


1  2  2 
W  =  r-  cot  d>  cosec  ip(2  +  sin  it)  cos  20  . 
0 


For  the  case  when  the  patch  is  centred  at  6  =  45°,  $  =  90°  these  displacements 
may  be  written  as 


<  .2  .4 

=  cos  <m-  (4-2  sin  <j>  +  sin  ij>)  cos  0  cos  20 


+  (2  +  sin2ij>)(2  sin  0  sin  20  +  sin  <f>  cos  0  cos  20)  sin^4>  , 
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<p  | -(4  -  2  sin^<j>  +  sin%)  sin  0  cos  26  1 

-  (2  +  sin2<}>)(2  cos  0  sin  !0  -  sin  <f>  sin  0  cos  20)j^6  sin3<(>  ,V  (4—  1 3) 


(1  +  cos  <j>)  cos  20  , 


2  sin 


(  .2  4 

(o  ■  (4  -  2  sin  <j>  +  sin  <$)  sin  cos  20 

<f>  \ 

+  3(cos2<j>  +  1)  cos  20  sin^cf  , 


2  *  *4 

0)  =  cos  4>  (2  +  sin  <j>)  ( 1  -  sin  <J>)  sin  20/3R  sin  <\ 

0 


The  quantities  U  ,  U  ,  U  are  prescribed  on  double  precision  SPC  cards 
x  y  z 

on  CPO.  co  and  go  are  resolved  into  co  ,  to  and  go  ,  see  equation  (4-7) , 

0  9  x  y  z 

and  prescribed  on  double  precision  MPC  cards  on  CPI.  Once  again  it  is  necessary 
to  introduce  scalar  points,  as  described  in  subroutine  KTHEPHI,  and  these  dis¬ 
placements  are  fixed  on  single  precision  SPC  cards  output  to  CP2. 

4.5  Subroutine  ETHETHE 

This  subroutine  is  called  when  ICASE  =  4.  It  prescribes  the  boundary 
displacements  for  the  nearly  constant  strain  state 


where 


—  constant  , 


(4-14) 


e„.  —  0 


The  displacements  are 


~  J^3  +  (2  +  cos  $)  cos  <{>  |  tan2  ^ 

-  1 3  -  (2  -  cos  <J»)  cos  <)>  j  cot2  -|j  cosec  <J>  cos  20  , 


2  2 
-  —  cot  4>  cosec  <)>  cos  26  , 


(4-15) 


2  3.2 

U  »  -  rr  cosec  ip  sm  6  * 

0  J 

1  2 

W  =  -j  cosec  4>  cos  20 


12 


When  the  patch  is  centred  at  9  ■  0°, 


90 


these  may  be  written  as 


2  3 

**  (3  sin  <j>  cos  0  cos  20-2  cos  38) /3  sin  <j>  , 

2  3 

Uy  *  (3  sin  $  sin  0  cos  20-2  sin  30) /3  sin  ip  , 

2 

U  *  cos  ij >  cos  20/sin  $  , 
z 

3 

0)  *  -  4  cos  ip  cos  20/3R  sin  ip  , 


>  (4-16) 


w0  “  0  ‘ 


J 


The  quantities  U^,  U  ,  are  prescribed 

on  CPO.  (o„  and  w,  are  resolved  into  u  ,  ui 
6  <p  x y 

and  are  prescribed  on  double  precision  MPC  cards 
introduced,  as  described  in  subroutine  KTHEPH1, 
single  precision  SPC  cards  output  to  CP2. 


on  double  precision  SPC  cards 

and  u  ,  as  in  equation  (4-7) 
z 

on  CPI.  Scalar  points  are 
and  their  displacements  £ixed  on 


5  NASTRAN  DATA  GENERATION  PROGRAMS 

The  programs  in  this  section  are  to  convert  free  format  files  containing 
either  the  coordinates  of  a  point  or  the  node  numbers  around  an  element  into 
output  that  can  be  directly  incorporated  into  a  NASTRAN  bulk  data  deck. 

5. 1  Program  GENCQUAD4 

This  program  gives  the  element  numbering  around  a  QUAD4  quadrilateral 
plate  element  (which  is  a  warped  isoparametric  plate  membrane-bending  element). 
The  output  is  in  the  form  of  NASTRAN  CQUAD4  cards. 


The  input  is 


N  N 1  N2  N3  N4 


where  N  is  the  element  number. 

Nl,  N2,  N3,  N4  are  the  node  numbers  around  element  N  .  This  is  repeated 
for  each  element  until  the  end  of  the  input  file  is  reached. 

The  output  in  fixed  format  for  a  constant  thickness,  homogeneous  isotropic 
material  is  in  the  form. 

CQUAD4  N  1  Nl  N2  N3  N4  0  . 

The  third  field  refers  to  the  property  card  and  the  eighth  field,  to  the  material 
orientation. 
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5.2  Program  GENCTRIA3 

This  program  outputs  CTRIA3  cards  that  can  be  included  directly  into  a 
NASTRAN  bulk  data  deck.  CTR1A3  cards  are  used  to  define  a  triangular  plate 
element  (compatible  with  the  QUAD4  quadrilateral  element). 

The  input,  all  integers  in  free  format,  is 

N  N1  N2  N3 


where  N  is  the  element  number, 

Nl,  N2,  N3  are  the  node  numbers  around  node  N  . 

This  last  line  is  repeated  until  the  end  of  the  input  file  is  reached. 

The  output  is 

CTRIA3  N  1  Nl  N2  N3  0. 

The  value  of  1  in  the  third  field  refers  to  the  NASTRAN  property  card  and  0. 
in  the  last  field  gives  the  orientation  of  the  material,  ie  for  a  homogeneous 
isotropic  material. 

5.3  Program  GEN GRID 

This  program  takes  as  input  the  node  number  and  corresponding  x,  y  and 
z  coordinates  of  a  mesh  and  outputs  them  in  the  form  of  GRID  cards  for  including 
in  a  NASTRAN  bulk  data  deck.  In  addition,  there  is  the  option  of  specifying 
freedoms  that  are  constrained  to  zero,  ie  permanent  single  point  constraints. 

The  input  data  is  read  in  in  free  format  and  is 

NNODE 

N.  X.  Y.  Z. 

1111 

Nj  IFREEDOM 

where  NNODE  is  the  total  number  of  nodes , 

X.,  Y.,  Z.  are  the  x,  y  and  z  coordinates  at  node  i  , 

ill 

(i  -  1,  2,  ....  NNODE)  , 

N.  is  the  node  at  which  the  freedom  is  being  constrained, 
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IFREEDOM  is  a  number  consisting  of  any  of  the  integers  1  6  with  no 

embedded  blanks,  the  integer  contains  1  if  Uxj  ■  0  , 

2  if  Uy.  -  0  , 


3 

4 

5 

6 

is  repeated  for  each  node  to  be  constrained  in  this  way. 
The  output  is  either  of  the  form 


if  Uy. 

3 

if  Uz.  ■  0  , 

J 

if  o>Xj  ■  0  , 

if  wyj  ■  0  , 

if  uz.  ■  0  , 

J 


or 


GRID 

I 

blank 

X. 

l 

Y. 

l 

Z. 

l 

GRID 

I 

blank 

X. 

l 

Y. 

l 

z. 

1 

NB  If 


c. ,  y .  or  z.  «  -1 
l  i  i 


the  output  field  width  is  exceeded  and  the  number  is  flagged  by  an  asterisk. 

The  output  file  must  then  be  edited  to  reduce  the  number  of  eight  characters 
and  any  equal  positive  numbers  reduced  similarly  to  maintain  symmetry. 

6  PROGRAM  CHANGEDATA 

This  program  converts  SEMILOOF  data  to  'facet'  data  by  altering  the 
position  of  the  midside  nodes  to  be  halfway  along  the  straight  line  joining  the 
two  vertex  nodes.  It  takes  as  its  input  a  SEMILOOF  data  file  and  outputs  the 
node  number  and  new  coordinates  to  a  card  punch  file  that  can  then  be  edited  back 
into  the  original  file.  If  there  are  more  than  200  nodes  in  the  data  file,  the 
program  has  to  be  altered  to  increase  the  size  of  arrays  X,  Y  and  Z  . 

7  CONCLUSIONS 

The  computer  programs  presented  here  are  a  faster  and  more  accurate  method 
of  data  generation  than  manual  preparation. 

The  programs  have  been  used  extensively  to  create  all  of  the  data  required 

.2 

for  the  problems  considered  by  Morley  and  Morris  . 
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SEMILOOF  PROGRAM 


The  program  listed  below  generates  the  element  stiffness  matrix  and 
associated  matrices  for  the  SEMILOOF  finite  element.  The  master  segment 
reads  in  the  data  and  then,  for  each  element  in  turn,  calls  subroutine  ELSTIF, 
which  calls  other  subroutines  to  calculate  the  element  stiffness  matrix. 
Immediately  before  control  is  returned  from  subroutine  ELSTIF,  the  matrices 
needed  to  calculate  the  stresses  are  written  to  channel  8.  On  returning  to  the 
master  segment,  the  element  stiffness  matrix  is  written  to  channel  7  and  the 
matrices  containing  pressure  and  gravity  loads  are  written  to  channel  9.  The 
next  line  of  element  connectors  are  then  read  and  the  matrix  calculation  process 
repeated. 


non  <->  o  n  n  o  o  non  non 
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SHORTLIST 
PROGRAM(LOOF) 

INPUT5=CR0 
OUTPUT  6*LP0 

0UTPUT7=E07/UN FORMATTED <DAF?) 

0UTPUT8  =  E  DA /UN  FORMATTED  (DA  FT.) 

0UTPUT9  =  ED*?/UNC0RMATTED(DAP°' 

TRACE  0 
END 

MASTER 

C****  PROGRAM  TO  GENERATE  THE  ELEMENT  STIFFNESS  MATRICES  AND  ASSOCIATED 
O***  MATRICES  FOR  THE  SEMILOOF  SHELL  FINITE  ELEMENT 
COMMON/ FR0n/ELST<528> ,ELR(37.2) , LN0DS<8> . LNODZ 
COMMON /SHELL/ ARE A.  ELXYZT(9,4) , F RAM ( 3 , 3) , PO I  NT ( 3) , 

1  SIDE,  THIK.  WSHEL(13,45)  XITA(2),  XYZ(241,3> 
COMMON/ERNFST/YMODT (10) ,  POIST(IO),  THI KT ( 1 0) »  DENSTT(IO), 

1  PRESST  (10)  ♦  I.MAT  ( 241 )  ,  STreSM(24 , 32) 

DO  800  1=1,13 
DO  800  J=1 ,45 
USHELd  ,  J)«0.0 
800  CONTINUE 


****  READ  probifm  parameters. 


4  READ ( 5  »  500>  NELZ 
500  FORMAT (10) 

WRITE(6,601 )NELZ 

601  FORMAT  ( /7H  NELZ  =,I4,<>H  ELEMENTS) 

****  READ  MATERIAL  PROPERTIES  AN,)  NODAL  COORDINATES. 

MATZ=1 

DO  6  MAT  =  1 ,MATZ 

****  YMODT  CONTAINS  youngs  MODULUS 
****  POIST  CONTAINS  POISSONS  R AT  1 0 
****  THIKT  CONTAINS  SHELL  THICKNESS 
****  DENSTT  CONTAINS  MATERIAL  DENSITY 
****  PRESST  CONTAINS  PRESSURE  LOADING 

READ(5,502)  VMOdT(MaT),  POIST(MAT).  THIKT(MAT),  DENSTT(MAT), 

1  PRESST(MAT) 

502  FORMAT(EO.O,4FO.O) 

6  CONTINUE 

****  NZ  CONTAINS  THE  TOTAL  NUMBER  OF  VETEX  AMD  MIDSIDE  NODES 
READ(5,591)  NZ 
591  FORMAT ( 10) 

****  FOR  EACH  NODE  READ  IN  NODE  NUMBER  AND  X,Y,Z  COORDINATES 
READ(5.504)  { N  ,  ( XYZ ( N , I ) ,  I  =  1,3),  N  =  1,NZ) 

504  FORMAT(IO,3FO.O) 

WR I TE ( 6 , 604 )  ( N ,  (XYZ(N,I>  I  =  1,3).  N  =  1.NZ) 

604  F0RMAT(/3X,1 HN,6X,1 HX,9X,1 HV,9X.1HZ/(1X. I3,3F10.6)> 

WR  I  TE ( 6 , 606) 

606  FORMAT  (  / 1  5  H  ELEMENT  NUMBER  ,  :,X  ,  20HNODE  NUMBERS  a  LNODS) 

****  FOR  EACH  NODE , READ  I N  ELEMENT  CONNECTIONS  AND  CALL  SEMILOOF  ROUTINES 

DO  7  NEL  =  1  .  NFt.Z 
LMAT(NED=1 
N  ■  NEL 

READ ( 5  »  5?2 )  (  LNODS  ( I )  ,  I  =  ^,8) 

597  FORMAT (810) 

WR  I  TE ( 6 ,608)  N , ( LNODS (I), 1*1, 8) 
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608  F0RHAT(2I11 ,915, (122,915)) 

C****  LNODZ  «  6  FOR  A  TRIANGULAR  ELEMENT  OR  s  8  FOR  A  QUADRILATERAL  ELEMENT 
LN0DZ*8 

I F  <  LNODS ( 8 ) . EQ.O)LNOnZ»6 
C 

C****  CALL  SUBROUTINE  TO  SET  UP  ELEMENT  STIFFNESS  MATRIX  AND  ASSOCIATED 
O***  MATRICES 

CALL  ELSTIF(NEL) 

C****  WRITE  LNODZ  CONTAINING  NUMBER  OF  NODES  AROUND  THE  ELEMENT  AND  ELST 
O***  CONTAINING  THE  ELEMENT  STIFFNESS  MATRIX  TO  CHANNEL  7 
WR I TE  <  7) LNODZ 
WRITE  (7)EL<!T 

C****  SEE  SUBROUTINE  ElSTiF  FOR  PRINT  OUT  OF  LNODZ , STRESM , PT , FR 
C****  CONTAINING  STRESS  INFORMATION  TO  CHANNEL  8 
C****  WRITE  LNOD7  AND  ELR  CONTAINING  PRESSURE  AND  GRAVITY  LOADS 
C****  TO  CHANNEL  9 
WRITEC9) LNODZ 
WRITE (9)<ELR(I ,1 ) ,1=1 ,32) 

7  CONTINUE 
STOP 
END 

BLOCK  DATA 


«***  TO  INITIALIZE  COEFFICIENTS  FOR  CORNER-MIDSIDE  AND  LOOF  VERSIONS 


****  OF  QUADRATIC  TRIANGLE  AND  QUAPR I LLATERAL  FOR  SUBROUTINE  SFR 


DIMENSION  COFFA< 

common/coef/coef 

EQUIVALENCF  (COF 
DATA  COEFA/  1..- 
1  -1 . ,  0. ,  2. .  0. 

2  2*,  0.,  0.,  4 . 

3  -6.04145lP84.-6 

4  -0.244016034,  0 

5  -4.928203230,-1 
4  -1.464101615,  5 

7  -2.577350269,-1 

8  -1.732050808,-0 

9  -4.464101615,  4 

1  1.577350269,  p, 

2  -1 . ,6. ,6. , -6 • » - 

3  .5,0. ,-.5, -.5,0 

4  .5,.5,0.,0.,0.f 

5  .5,0. ,.5, -.5,0. 

6  .5»-,5»0.,0..0 
DATA  COEFB/  0-0 •  j 

1  0.216506351,  o 

2  -0.000000000,-0 

3  0.281250000,  0 

0.375000000,  0 
-0. 375000000, -0 


,,  0. 
0.  , 


0. 


166).  COE  F  B ( 81 > 

(2  47) 

F ( 1 ) ,  COEFA < 1 ) ) ,  (COEF(167> ,COEFB(1)> 

3., -3.,  2.,  4.,  2.,  0.,  4.,  0 . , -4 . , -4 . ,  0 
,  0.,  0.,  4.,  0.,  0,,  0 . , “1 . , 

,  0.910683603,  1.577350269, 

2.464101615,  8.928203230,  1.732050808, 
2.041451884,  4.196152423,-4.464101615, 
0.333333333,-1.422649731,-2.577350269, 

5.464101615,  1.732050808,  0.333333333, 
5.000000000,-1.464101615, 
0.422649731,-4.928203230, 
0.910683602,-6.041451884. 
6.196152422,-1.732050807 

25, -.25, -.25,0. , 

25, -.25, .25, .25, -.25,0.  , 


2 

,  0. ,  0.,  0 
,  0 . , -4 . , -4 
.196152423, 
.422649731 , 
.732050808, 
.000000000, 
.422649731 , 
.244016936, 
.196152423, 
.928203230, 
6, 6. ,0. , — 
.,0..0.,.5,0 


5.464101615, 
2.041451884, 
1.732050807, 

2.464101615, ' 


25, 

r  0  •  ,  0  •  , 

.25, 

.25, 

o . » ■ 

-.25,0. 

,0.  , 

.25, 

• .  ,  • 

-.25,0. 

1 0 1 1 

.25, 

0 . ,  ■ 

-.25,0. 

9  0  •  9 

.25, 

0.  , 

■  1  .  ,  0  .  , 

o. 

1  .  ,0 

-0.216506351  ,  (1 
0.649519053,  0 
-0.093750000,  0 
-0.324759526,  ,i 
-0.216506151  ,  1 

2  -0.000000000,-0 

3  -0.093750900,  0 

4  -0.375000000,  .) 

5  .375, -.64051905 
END 


— .5,— .5,0.  ,  —  .25,0.  ,0.  ,  .25,  .25,  .25,  .25,  .25,0. 

,0- ,0. , — . 5 , 0 . ,— .25,0. ,0* , .25,— .25, .25,— .25, .25,0, 

,— 1.,0.,0.,1./ 
0000000,  0.2165)6351 ,-0.375000000,-0.093750000. 
.281250000,-0.649519053,  0 . 375000000 , -0 . 324759526 , 
.216506351  ,-0.375000000, -0.093750000, -0.21  6506351  , 
.649519053,  0.375000000,  0.324759526,  0.000000000, 
.216506351 ,  0.281250000,-0.216506351,-0.093750000. 
.649519053,-0.324759526,  O.oOOOOOOOO,  0.375000000, 
.281250000,  0.216506351 ,-0.093750000,-0.375000000, 
.324759526,-0.000000000,-0.216506351 ,  0.375000000, 
.216506351,  0.281250000,  0.649519053,-0.375000000, 
.000000000,  0.216506351,  0.375000000,-0.093750000 
.281  250000,-0. 649519053, -0.3750(10000,  0.324759  5  26, 
.375000000,-0.216506351 ,  0.281250000,-0.216506351 , 
.375000000,  0.649519053,-0.324759526,-0.000000000, 
.216506351,  0,281250000,  0.2165)6351,-0.093750000, 
3, .324759526.1 . ,0.,O.,-.75,0.,-.75,0.,0.,0./ 
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SUBROUTINE  ELSTlF(NEl) 

C 

C****  DEMONSTRATION  PROGRAM  FOR  ELEMENT  STIFFNESS  ETC. 

REAL  PT ( 3#  4) >  FR ( 3» 1 2) 

C 

DIMENSION  B(6<32),  BV(192>,  DM0D(6,6),  DMODV(36),  XGAUS(4,4), 
1  GVECT (6) *  GRAVY(3) 

COMMON/ FRON/ELST (528) ,ELR( 32,2) .LNODS (8) ,LNODZ 
COMMON /SHE  l  L/AREA  f  E  I.XYZT  (9.4)*  FRAM(3,3),  POINT(3), 

1  SIDE,  THIK*  W$HELM3,45),  X1TA(2>,  XYZ(241,3) 
COMMON/ERNEST/YMODT(10) ,  POlST(IO),  THIKT(IO),  DENSTT(IO), 

1  PRESST <1 0) «  lMAT  <241 ) ,  STreSM(24, 32) 

EQUIVALENCE  (GVECT  <1 ),  R)  ,  (GVECT (2) , S > ,  (GVECT ( 3) #  T) » 

1  (GVECT  <4) »U) ,  <GVECT<5> fV>,  (GVECT (6) ,U) , 

2  (B(1»1),BV(1)>#  ( DHOD ( 1 » 1 ) #  DMODV ( 1 ) ) 

****  initialise  arrays. 

data  XGAUS/0..  4*. 5,  0.,  2*. 5,  -.577350269,  2* . 577350269 , 

1  3*-. 577350269,  2* . 577350269/ 

LVABZ  «  4*LN0DZ 
IZ*(LVABZ*(LVABZ*1) ) /2 
DO  2  I  *  1 , IZ 

2  ELST(I)  «  0.0 

DO  3  I  *  1. LVABZ 

3  ELR(I,1)  «  0.0 
DO  4  N  «  1.LN0DZ 
DO  4  I  *  1.3 

4  ELXYZT(N,I)  ■  XYZ(LNODS(N),I> 

C 

C****  DEFINE  PROBLEM  CONSTANTS. 

C 

MAT  *  LMAT(NEL) 

YMOD  »  YMODT<MAT) 

POIS  *  POIST(MAT) 

OENSTY  *  DFNSTT(MAT) 

PRESS  *  PRFSST(MAT) 

THIK  •  THIXT(MAT) 

DO  6  N  *  1 . LNODZ 


6  ELXYZT (N 

f  4) 

9 

GRAVY ( 1 ) 

9 

0.0 

GRAVY (2) 

9 

o.o 

GRAVY (3) 

9 

1.0 

Appendix 
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C 

C****  CALL  HALOOF  AND  CREATE  MATRIX  B  AT  THE  INTEGRATING  POINTS. 

C 

IGAUSZ  >  LNODZ/2 
DO  26  IGAUS  =  1 , IGAUSZ 
XITA(1)  =  XGAUS(IGAUS,LNODZ-5) 

XITAC2)  *  XGAUS(IGAUS,LNODZ-4) 

CALL  HALOOF(NEL) 

DO  8  N  s  1 , LVABZ 
B ( 1 ,  N  )  i  WSHFLU.N) 

B(2*N)  a  WSHEL(7,N> 

B(3#N)  a  WSHEL(5.N)  +  WSHEL'6#N) 

B  ( 4  #  N  >  a  WSHELdO.N) 

B(5,N)  a  WSHEL(12,N) 

8  B  <6  *  N )  a  2.0*USHEL(11«N) 

C 

c****  create  modulus  matrix  for  calculating  the  stiffnesses. 

c 

DO  10  I  *  1,36 
10  DMODV(I)  ■  0.0 

GASH  a  YM0D*AREA*THIK/(1 . O-pO I S*P0l S > 

DMODV(I)  =  GASH 
DMODV(Z)  a  GASH*POIS 
DMODV <7)  a  GASH*POIS 
DMODV (8)  a  GASH 

DMODV (15)  a  0.5*GASH*(1 .O-POIS) 

GASH  «  GASH*THI«*THIK/12.0 

DMODV ( 22 )  a  GASH 

DMODV  < 23 )  a  GASH*P0IS 

DMODV (28)  a  GASH*P0IS 

DMODV ( 29 )  a  GASH 

DMODV ( 36 )  a  0.5*GASH*(1 .O-POIS) 

C 

C****  GET  DB,  A  COLUMN  AT  A  TIME,  AND  COMPUTE  STIFFNESSES. 

C 

NSTIF  «  0 

DO  24  NROW  a  1 ,  LVAB7 

DO  14  K  »  1<6 

GASH  a  0.0 
LDEL  a  6MNR0W-K) 

LA  »  6*K-5 
LZ  »  $*k 

DO  12  L  »  LA.LZ 
11*L*LDEL 

12  GASH«GASH+BV(L1) * DMODV ( L) 

14  GVECT(K)  >  GASH 
IDEL  *  -5 

DO  16  KOL  a  1,NR0W 
IDEL  ■  IDEl  ♦  6 
NSTIF1aNSTlF*K0L 

16  EIST(NSTIF1 )aELST(NSTIFl >+R*BV(IDEL)+S*BV(IDEL+1) 

1  ♦  T*BV  <  t  DEL  +  2)  ♦  U*BV ( I DEL+3)  ♦  V*BV(IDEL+4)  ♦  W*BV(IDEL*5) 
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C****  accumulate  pressure  and  gravity  LOADS. 

C 

GASH  *  ELRfNROw.l) 

DO  18  I  *  1,3 
GRAVY(I)  *  POINT(I) 

GRAVY (3)  *  0.0 

FACT  *  AREA*THIK*DENSTY*GRAVY(I> 

DO  18  J  =  1,3 

18  GASH  *  GASH  +  F ACT*FR AM < I , J ) *WSHE L ( J , NROW) 

ELR(NR0W,1)  *  GASH  +  PRESS*AREA*WSHEL(3,NROW) 

c 

C****  COMPUTE  AND  STORE  STRESS  MATRIX  AT  GAUSS  POINTS. 

C 

DO  22  I  *  1,3 
J  *6*  ( I GAUS-1 ) ♦ I 

STRESM(J,N»OW)  =  GVECT < I) / < AR EA*TH I K > 

22  STRESM(J  +  3, NROW)  *  GVECT ( 1  +  3 ) *6 . 0/ ( AR EA*TH I K*TH I K) 

24  NSTIF  =  NSTIF  ♦  NROW 
DO  7  1*1,3 

PT(I ,  IGAUS)=pOlNT(I) 

DO  9  J=1,3 
NP*3*(lGAUS-1 )+J 
9  FR(I,NP)=FRAM(I, J) 

7  CONTINUE 
2*  CONTINUE 

C****  WRITE  LNODZ, STRESM.pT , FR  CONTAINING  STRESS  INFORMATION 
C****  TO  CHANNEL  8 
WRITE(8)LN0DZ 
WRITE(8)STRESM 
WRITE(8)PT 
WRlTE(8)FR 

return 

END 

SUBROUTINE  HALOOF(NEl) 

C 

C****  TO  CREATE  SHAPE  FUNCTION  ARRAY  WSHEL,  FOR  SEMILOOF  SHELL  ELEMENT. 

C 

C****  WRITTEN  BY  BRUCE  IRONS,  JULY  1972, WASHINGTON  D.C. 

C 

DIMENSION  AREAV(3>,  FRAME (3,3),  GENSID<6,4), 

1  SHEARdl  ,43) ,  S I GT  ( 3)  ,  SWOP  ( 6) ,  THIKDD<3,3>,  TRANS(2,2), 

7  VLOOF (3,36) .  WCOrn (10,3) ,  WLOOF(10,3),  XGAUS(4,4).  XILOOF(9,4), 
3  XLOCA L ( 2) ,  XYZDD ( 3,3),  XVZPRE(8,4) 

COMMON /F RON /ELST (528) ,ELR(32 ,2) , LNOdS  <8> , LNODZ 
COMMON/SHEt L/AREA,  ELXYZT<9-4>,  FRAM(3,3),  POINT(3), 

1  SIDE,  THIK,  WSHELM3,45),  XITA(2),  XYZ<241,3) 

EQUIVALENCE  < T1 1 . TRANS ( 1 , 1 )> ,  <T1 2 , TRANS ( 1 , 2 ) ) , 

1  <T21,TRANS<2,‘i>>,  (T22, TRANS  (2,2) ) 

DATA  GENSlD/1 . ,  -1.,  0.,  3*-. 5,  0.,  1.,  -1.,  4*1.,  0.,  -1., 

1  4*0.,  1.,  0.,  -1.,  2*1. /XILOOF/. 211324866,  2* . 788675 1 34 , 

2  .211324866,  2*0.,  .3333333333,  4*0.,  .211324866,  2* . 7886751 34 , 

3  .211324866,  .3333333333,  2*0.,  -.577350269,  .577350269,  2*1., 

4  .577350269,  -.577350269,  2*-1 . ,  0.,  2*-1 . »  -.577350269, 

5  .577350269,  2*1.,  .577350269,  -.577350269,  0.  /, 

6  XGAUS/0 . ,  4*. 5,  0.,  2*. 5,  -.577350269,  2* . 577350269 , 

7  3*-. 577350264.  2* . 577350269/ ,  XYZPRE/32*0.0/ ,  NOZPRE/0/ 
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C 

C**++  GENERATE  NSTAGE  TO  DEFINE  PATH  THROUGH  HALOOF. 

C 

N  ERROR  =  1 

IFUNODZ.NE.A.  AND  .LNODZ. NE. 8)  GO  TO  90 
NSTAGE  *  4 

IF(LNODZ.NF.NOZPRE)  NSTAGE  =  2 
NOZPRE  >  LNODZ 
DO  2  LNOD  =  1  . 1 NODZ 
DO  2  NX  =  1  ,4 

IF(ELXYZT(LNOD.NX) .NE.XYZPRE(LNOD,NX))  NSTAGE  *  2 

2  CONTINUE 

IFCNSTAGE.FQ.4)  GO  JO  18 
C 

C****  INITIALIZATION  FOR  NEW  ELEMENT.  NSTAGE  a  1.  FIND  CENTRE  COORDINATES 
C 

LIMZ  *  <3* | NODZ) /2  -  1 
LVABZ  *  4* L NODZ 
LNODZA  >  LNODZ  ♦  1 
LVABZA  »  LVABZ  +  1 
LVABZZ  «  LVABZ  ♦  LIMZ 
DO  3  L  a  LNODZA, LIMZ 
DO  3  J  *  1 , LVABZZ 

3  SHEAR  <  L ,  J  )  a  0.0 
DO  5  NX  =  1,4 
GASH  »  0.0 
LNODZWNODZ/2 

DO  4  KORNai f LN0DZ1 
DO  4  K  a  1 .2 

K1 *2*KOrN+k-2 

K2*216*K-408- LNODZ* ( 21 *K-41 ) 

4  GASH  *  GASH 

1  +8.0*ELXY7T^X1 ,NX)/FlOAT<K2> 

•5  ELXYZT(0,NX)  »  GASH 
C 

C***  +  DIAGNOSTICS  FOR  A  NEW  ELEMENT.  RELATE  COORDINATES  TO  CENTRE. 

C 

DO  10  I  »  i ,  LNODZ 
N  ERROR  a  2 

IF(ELXYZT(t,4).LE.0.O)  GO  TO  99 

IFCI.EO. LNODZ)  GO  TO  9 

JA  ■  I  +  1 

DO  8  J  >  JA,  LNODZ 

N  ERROR  a 

IFCLNODS(I) .FQ. LNODS(J))  GO  TO  99 
DO  7  K  a  1 .3 

IF(ELXVZT(I,IC>.NF.ELXYZT(J,<>>  GO  TO  8 

7  CONTINUE 

N  ERROR  a  4 
60  TO  99 

8  CONTINUE 

9  DO  10  NX  «  1.4 

£  I F (NX.NE.4)  E LXYZT ( I , NX)  ■  dLXYZT(I.NX)  -  ELXYZT(9,NX) 

10  XYZPRE  <  I ,  NX)  a  ELXYZT(I#NX) 


U 
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C 

C****  CREATE  SWOP  *  1.0  OR  -1.0.  TO  IN 

C****  CREATE  SWOP  »  1.0  OR  -1.0,  to  IMPLEMENT  SIGN  CHANGES  AT  LOOF  NODES. 
C 

C***«  ALSO  INTERPOLATE  TO  ESTIMAT-  NORMAL  THICKNESSES  AT  I.OOF  NODES. 

C 

VLOOF ( 1 <  LVABZA)  *  ELXYZT(9.4> 

DO  12  NSIDF  =  1,6 

1 2  SWOP (NS  IDE)  b  1 .0 
LAST  *  LNODZ  -  1 

DO  14  NEXT  s  1.  LNODZ.  2 
MID  =  LAST  ♦  1 

IF(LNODS(NEXT).LT.LNODS(LAST>)  SWOP (MID/2)  =  -1.0 
LAST1«4*LA$T-3 

VLOOF  (1  ,LAST1>*. 455341 801  *E'.XYZT( LAST. 4) 

1  ♦  . 666666667* E LXYZT (MI D , 4 )  -  . 1 22008468*E LXYZT (NEXT , 4) 

MI D1 =4*MI D-3 

VL00F(1 .MIDI)  =  -.122008468*ELXYZT(LAST,4) 

1  *  .666666667#E LXYZT (MID. 4)  +  . 45 5341 801 *E LXYZT (NEXT . 4) 

C 

C****  ALSO  CHECK  THAT  MIDSIDE  NODES  ARE  REASONABLY  CENTRAL. 

C 

GASH  *  0.0 

GISH  *  0.0 

GUSH  =  0.0 

DO  13  I  *  1,3 

ELMID  «  ELXYZT(MID, I) 

GASH  «  GASH  ♦  (ELXYZT(NEXT.I)-ELMID>**2 

GISH  »  GISH  +  (ELXYZT(LAST,T)-ELMID>**2 

13  GUSH  «  GUSH  «■  (ELXYZT(LA$T, I >♦* LXYZT ( NEXT , I ) -E LMI D-E LMI D> **2 

N  ERROR  =  5 

IF(ABS(GASH-GlSH). GT. 0.040*'. GAS H+GISH))  GO  TO  Q9 
N  ERROR  =  6 

C  IF(GUSH.GT.0.25*(GASH+GISH) )  GO  TO  99 

14  LAST  =  NEXT 
C 

C****  ORGANISE  LOOP  AROUND  NODES,  FOR  NSTAGE  b  2 
C 

C****  DO  76  NSTAGE  *  2,4  (IN  NSTAGE) 

15  NLOOF  *  0 

16  NLOOF  »  NLOOF  ♦  1 

C****  DO  67  NLOOF  s  1 ,  LNODZ+1  IF  NSTAGE  b  2, 

C****  OR  DO  67  NLOOF  a  1,  (3*LNOr>7)/2  IF  NSTAGE  »  3. 

DO  17  I  =  1,2 
LN0DZ2*LN0dZ+I-6 

I F (NSTAGE . FQ.2 .  OR  . NLOOF . LE . LNODZ) 

1  XLOCAL(I)sXlLOOF(NLOOF,LNODZ2) 


C 

C****  AND  ALSO  Around  integrating  POINTS  IF  NSTAGE  s  3. 

c 

NLOOF1*NLOOF-LNODZ 

LN0DZ2«LN0dZ*I-6 

IF(NSTAGE.EQ.3.  AND  . NLOOF. <',T.  LNODZ) 

1  X  LOCAL ( I ) eXGAUS (NlOOFI , LNODZ 2) 

17  CONTINUE 
GO  TO  23 
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****  OTHERWISE,  ORGANISE  SINGLE-SHOT  OPTION,  FOR  NSTAGE  *  4. 

***★  TEST  WHETHER  INPUT  POINT  IS  A  LOOF  NODE,  PLUS  OR  MINUS  0.0001. 

18  DO  19  I  =  1 ,2 

19  XLOCAL(I)  =  XITA(I) 

NLOOF*  lnodza 

DO  22  MAYBE  =  1  ,  LNOdZ 
DO  20  I  »  1  ,2 

IF(ABS<XLOrALO)-XlLOOF(MAYHE,LNODZ+I-6)>  .GT. 0.0001)  GO  TO  22 

20  CONTINUE 
NLOOF  ■  MAYBE 

22  CONTINUE 

****  CREATE  VALUES  AND  XI, ETA  DERIVATIVES  OF  X,Y,Z  IN  XYZDD ,  T  IN  THIKDD 

23  CONTINUE 

CALL  SFR (XLOCAL ,  WCORN,  WLOOF,  NSTAGE) 

K  a  0 

DO  27  1=  1,3 
DO  26  J  a  1,3 
GASH  >  o.O 
DO  24  L  *  1  , LNODZ 
L1-L+K 

24  GASHaGASH+wcORNUl  ,  I)*ELXYZT(L,  J) 
xyzdd ( J, i )  «  gash 
IF<NSTAGE.eq.2)  GO  TO  26 

GASH  a  o.O 

DO  25  L  a  1, LNODZA 

Li  * L*K 

L2*4*L-1 

25  GASH«GASH*WLOOF(Ll , I ) *VLOOF < J , L2> 

THIKDD(J.I)  a  GASH 

26  CONTINUE 

27  K  «  1 

****  create  vector  area  a  varea,  at  given  point  xi ,  eta. 

CALL  VECT0r<XYZDD(1  ,2) ,  XYZnD(1,3>.  AREAV(D) 

CALL  SCALAR (AREaV<1 ) ,  AREAV(I),  AREASQ) 

N  ERROR  *  7 

IF(AREASQ.EQ.O.O)  GO  TO  99 
AREA  ■  SQRT(AREASQ) 

****  NORMALISE  VECTOR  AREA  INTO  FRAME#  COL. 3,  AS  LOCAL  UNIT  NORMAL  Z. 
****  COLUMN  2  OF  FRAME  BECOMES  UNIT  Y  AROUND  EDGE. 

DO  30  I  a  i  ,3 

FRAME (1,3)  a  ARE AV ( I ) /AREA 
GASH  ■  0.0 
DO  29  J  a  1  ,2 
NLOOF1a(NLOOF*1)/2 
LN0DZ3*LN0dz*J-6 

29  GASHaGASH*GENSlp(NLOOF1 , LNOoZ3) *XYZDD < I  ,  J  +  1 ) 

30  FRAME (1,2)  a  GASH 
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C 

o***  normalise  Y,  AND  IMPLEMENT  swop  by  reversing  sign  of  y. 

c 

o***  PUT  APPROXIMATE  VECTOR  THICKNESS  ETC.  into  VLOOF,  FOR  NSTAGE  *  2 
C 

N  E R ROR  s  8 

CALL  SCALAR(FRAME<1 ,2) ,  FRAME (1.2).  SIDESQ) 

IF(SiDESQ.EQ.O.O)  GO  TO  09 
SIDE  *  SORT (SIDESQ) 

DO  31  I  ■  1.3 
NL00F4*(NL00F+1 ) /2 

FRAME(!.2)sFRAME(!’2>*SUOP<NLOOF4> /SIDE 

IF<NSTAGE.NE.2>  GO  TO  31 

NL00F5»4*Nl OOF-2 

VLOOF ( I , N  L00F5 ) =FRAME(I»2) 

NL00F6«4*Nt00F-1 
N L00F7*4*N I  OOF-3 

VLOOF (I,NL00F6) =FRAME(I #3)*VL00F<1 .NL00F7) 

NL00F8»4*NL00F 

VLOOF(I,NLOOF8)=FRAME(I,3> 

31  CONTINUE 
C 

C***»  AND  COLUMN  1  IS  UNIT  X,  THE  OUTWARD  POINTING  IN-PLANE  NORMAL. 

C 

CALL  VECT0r<FRAME<1,2),  FRAME(1,3).  FRAME<1,1>) 

C****  CHECK  THAT  NORMALS  ARE  REASONABLY  PARAPPEL,  WHILE  NSTAGE  =  2. 

C 

IF(NSTAGE.GT.2)  GO  TO  35 
IF(NLOOF.EO.I)  f,0  TO  67 
KZ  *  4*NL00F-4 
DO  32  K  *  4,  K2,  4 
NL00F9=4*NL00F 

CALL  VECTOR (VLOOF (1 ,NL00F9> .VLOOF (1 , K>  » POINT (1 > > 

CALL  SCALAR(P0lNT<1)  ,  POINT'D,  COSSQ) 

N  ERROR  =  9 

C  IFCCOSSQ.GT.0.75)  GO  TO  99 

32  CONTINUE 
C 

C****  PLACE  CONTRIBUTION  OF  CENTRAL  NODE  IN  VLOOF  (NSTAGE  *  2  ONLY) 

C 

C****  COMPLETE  LOOP  NlOOF  »  1  TO  tNODZ+1  FOR  NSTAGE  =  2. 

C 

IF(NLOOF.LE.LNOpZ)  GO  TO  67 
THIKC  »  VLOOF ( 1 . LVABZA) 

DO  33  I  *  1 ,3 
DO  33  J  »  1,2 
LVABZ1*LVABZ+J 

33  VL00F(I,LVABZ1)=FRAME(I, J)*TH1KC 
GO  TO  67 

C 

C****  CREATE  THE  2X2  JACOBIAN  MATRIX.  AND  INVERT  IT.  (NSTAGE  =  3  OR  4) 
C 

35  DO  36  J  *  1 ,2 
DO  36  I  *  1  ,2 

CALL  SCALAR(FRAME(1 ,  I) .  XYZDD (1.J+1),  TRANS(J.I)) 

36  CONTINUE 
GASH  ■  Til 

Til  *  T22/AREA 
T22  *  GASH/AREA 
T1 2  «  -T12/AREA 
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N  ERROR  a  10 
IF(THIK.LE.O.O)  60  TO  99 
C 

O***  FIrtD  THE  CHANGE  IN  LOCAL  X.Y  DERIVATIVES  ACROSS  THICKNESS  OF  SHELL. 
C 

DO  57  LNOD  a  1 , lNODZA 
IF<NSTAGE.NE.4>  GO  TO  51 
DO  50  I  ■  ?,3 
GASH  «  0.0 
DO  49  J  *  1,2 

49  GASH  ■  GASH  -  TH!KDD<J.I)*WCORN(LNOD»J+1> 

50  POINT!!)  ■  GASH 
C 

C***«  CREATE  WSHFL  »  SHAPE  FUNCTION  ARRAY,  DISPLACEMENT  TERMS  FIRST. 

C 

51  KORN  ■  (LN0D+D/2 
DO  54  K  ■  1,3 

KOL  ■  2*K0rn  ♦  3*LNOD  ♦  K  -  5 
I  F  <  LNOD . GT  ,  LNOdZ >  KOI.  ■  5*LN0DZ  ♦  2  ♦  K 
DO  53  N  «  1,3 
FACT  •  FRAME(K.N) 

WSHEL(N,KOl)  ■  UCORN(LNOD,1)*FACT 
IF(NSTAGE.EQ.4.  AND  .N.EQ.3)  FACT  ■  0.0 
DO  53  NO  ■  2,3 
N3BN+N+ND 

53  WSHE L ( N3 »KOL)»«CORN( LNOD, ND)*F ACT 
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00  54  N  *  1,2 
00  54  NO  *  2.3 

USHEL(N+7,K0L)  =  WSHEL(N*7 »  <0L) 

M3*N+N+MD 

1  -THIKDD(ND-1,1)*WSHEL(N3»<0L)/THIK 
N4*N+N*ND+6 

I  ?  (NSTAGE  . FQ . 4)  WSHEKN4.  K0 !.  >  =  (  POINT  (ND) * FRAME (K, N> 

1  ♦  THlKDD(3/ND)*WCORN(LNOo.N+1)*FRAME(K#3))/THIK 

54  CONTINUE 

****  INTRODUCE  ROTATION  TERMS  WITH  BENDING  ACTION  INTO  WSHEL. 

DO  57  L  =  1,2 

KOL  =  (  L-1 ) *4*  LNODZ  *  (2-L>*6*K0RN  +  LNOO 
IFILNOD.GT. LNODZ)  KOL  ■  5*LN0DZ  ♦  3  -  L 
00  56  N  =  1,2 
LN0Dl*4*LN0D+L-4 

CALL  SCALAR (VLOOF  <1 , LN0D1 ) ,  FRAME (1 .N> .FACT) 

WSHEL  (N  +  7 ,  tfOL )  =  F ACT*WLOOF  f  »  NOo  ,  1 ) /TH I K 
IFINSTAGE.NE.4)  GO  TO  56 
00  55  NO  *  2.3 
N4=N+N+ND+6 

55  W$HEL<N4.K0L)=FACT*WL00F<LN00,ND)/THIK 

56  CONTINUE 

00  57  NROW  =1,7 

57  USHEL(NROW,kOL1  =  0.0 

****  COMBINE  LAST  THREE  COLUMNS  uF  WSHEL  TO  CREATE  NORMAL  DEFLECTION. 

I F ( LNODZ . Eo . 6)  GO  TO  61 
IZ  =  3*NSTAGE  +  1 
DO  60  I  =  1  ,  IZ 
,  GASH  *  0.0 
DO  59  K  s  1,3 
LN0DZ2=4*LN0DZ*4 

59  GASHsGASH+WSHEU  I  , K+42)*VLOOF  (K,  LNODZ2) 

60  WSHEL( I ,43)  =  GASH 

61  I F (NSTAGE . FQ. 4)  GO  TO  86 

****  CREATE  ARRAY  SHEAR.  FOR  INTRODUCING  THE  CONSTRAINTS  (NSTAGE  «  3) 

IF(NL00F.GT. LNODZ)  GO  TO  63 
DO  62  I  *  1 , LVABZZ 
SHEAR(NL00F,I)  =  WSHF.L<9  .  I ) 

SHEAR  (11.1)  =  SHEAR(H.I)  ♦  WSHE L (8 , I ) *S I DE*TH I K*SWOP ( ( NLOOF+1 > / 2) 

62  CONTINUE 
GO  TO  67 

63  DO  66  KOL  =  1 , LVABZ2 
00  66  NXY  =  1,2 
LN0DZ3*LN0r>Z*NXY 
GASH*SHEAR  < LNOD73.KOL) 

00  65  MXY  =  1,2 
LNOOZ4«4*LNOTZ+NXY 

CALL  SCALAR ( FRAME (1 ,MXY) ,VLO0f (1 , LN0DZ4) , FACT) 

65  GASH  *  GASH  ♦  WSHEL (MXY+7 »KOL)*AREA*THlK*FACT 
LNO0Z3«LNOnZ-*’NXY 

66  SHEAR  (LN0D73,K0i.)*GaSH 
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****  COMPLETE  LOOP  AROUND  LOOP  NODES  ETC.  TO  CREATA  VLOOF  OR  SHEAR. 

67  I F (NLOOF , LE . INOdZ.  OR. 

1  ( NSTAGE . EQ . 3 .  AND  .  NLOOF . LT . <3*LN0DZ) /2) >  GO  TO  16 

IF(NSTAGE.NE.2)  GO  TO  76 

****  CREATE  PLUS-MINUS  SUM  OF  THIKNESS  VECTORS  AT  LOOF  NODES  (NSTAGE>2) 

DO  70  I  »  1,3 

GASH  *  0.0 

DO  68  N  *  3  ,  LVABZ  »  4 

68  GASH  *  -GASH  ♦  VLOOfUiN) 

SIGT(I)  *  GASH 

****  AND  THE  3X3  MATRIX  ASSOCIATED  WITH  IT,  STORED  IN  XYZDD. 

DO  70  J  «  1,3 
GASH  »  0.0 

IFd.EQ.J)  GASH  =  FLOAT(LNODZ) 

DO  69  N  *  2, LVABZ, 4 
60  GASH  ■  GASH  -VLOOF < I , N) *VLOOF < J , N) 

70  XYZDD ( I , J )  «  GASH 

****  GET  THE  ADJUGATF.  OF  THIS  3X3  SYMMETRIC  POSITIVE  DEFINITE  MATRIX. 

K  «  3 

DO  71  I  *  1,3 
XI *6- I -K 

CALL  VECTOR (XYZDD <1 , 1 ), XYZDD <1 .<1 > , FRAME <  1  , K) > 

71  K  *  I 

CALL  SCALAP(XYZDD<1,1),  FRAME(1,1),  DETERM) 

DO  73  I  «  1,3 

CALL  SCALAR(FRAME<1 , T) ,  SIGT<1>,  PROD) 

73  POINT(I)  *  PROd/DETERM 

****  CORRECT  VECTOR  THICKNESS  IN  VLOOF. 

FACT  ■  1.0 

DO  75  N  ■  2  ,  LVABZ , 4 

FACT  ■  -FACT 

CALL  SCALAR (POINT (1 ) ,  VL00F(1,N>,  PROD) 

DO  74  l  *  1,3 

74  VL00F(I,N*1)  *  VLOOF  < I , N+1 )  -  FACT*(POlNT(I) -PROD* VLOOF ( I »N>  > 

****  CREATE  DIFFERENTIAL  DISPLACEMENT  VECTORS  TO  DEFINE  ROTATIONS. 

****  THIS  COMPLFTES  WORK  FOR  NSTAGE  ■  2 

TFIRST  ■  VLOOF(1,N-1> 

CALL  VECT0R<VL00F<1 ,N) ,  VLOOF ( 1 , N+1 ) ,  VLOOF ( 1 , N-1 ) ) 

DO  75  I  *  1,3 

75  VLOOF(I,N)  ■  VL00F(I,N)*TFIHST 
NSTAGE  ■  3 

GO  TO  15 
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C 

C*  +  *  + 

c 

c**** 

c 

76 


80 

C 

C***  + 

c 

c  **#* 
c 


"AS  "SEN  CUTS,  „  „lOW  l00p  f0B 
*  CH00SE  "'VOT  FOB  REOUCINS  ARRAY  SHrar 

.  continue  A*-  a’"‘  00  ,ow  ""“cnanee. 

??  83  UH  =  l.UMZ 
KP  3  tVASZ  ♦  L IM 
PIVOT  »  0.0 

JJ.J?.1  *  UMaLIMZ 

18 1 G  »(l  V0T> *6T.ABS<SNEARCL,KP>)>  <jo  TO  79 

change*  V’LVAB22 

S««fir.rBEARU8I6'lt> 

*«BARtllH.;f  -“cSiJM/pfJof* 

«8»«CE  ARRAY  Shear  T0  CREATE  * 

TH.S  COMPLETES  UORR  EOR  NST.LSE  .  3. 


c 

c**+* 


87 

88 
C 

C  *  *  *  * 

c 


91 

9? 

C 

C*  **  * 

c 


00  82  NROW  .  1 , L IH2 
P4CT  «  SHEAR CNROW/ Kp) 

IFINROW.EQ  1  yM  no  Cjl__ 

00  81  KOL  ;  1,LVABZ2 MCT,EQ  0,0)  60  T°  82 

81  SHEAR(NROW,<0L)  - 

82  CONTINUE  SHEAR(NROM,  Wl)  -  FACT*SHEAR  < ,  rn,  . 

8?  CONTINUE  "tAR(LiM,fC0U 

NST AGE  *  4 
GO  TO  18 

"  ““  A""AV  SHEA“  T0  CONSTRAIN  USHEL  AT  THE  G,VE»  ,» 

IS  00  as  ,  .  ,,lVA82  THE  G,VE»  Roint  xi.eta. 

DO  88  J  =  if13 
GASH.  WSHEL(J,D 
DO  87  K  a  1 ,  I M 7 
CVA8ZK«(C+LVAB  2 

8  Ss«uJ?;;“r!i<^CvA8«,.SHEA«(,.I) 

*  IHPIEHENT  swop  t0  exchange  tu0  norhal  ^ 

DO  92  N  *  8,LVAft7  » 

60  T0 

HSHEUJ.n-u  . 

CONTINUE  change 

ASSEMBLE  Uxz.  UYZ.  vxz  W7  TA 

’  VYZ  to  CREATE  WXX,  WXY,  WYY 

’?« /  i.ivabz 

VSHEU11 'U!  "  -f;ELMO,N) 

WSHEUIZe'n)  l  l2sHEU1S3EN)ni*NUWSHEl<12'N)) 

CONTINUE640'  'WS«*U;.N>.  J  s  1f12) 
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C 

C**«*  PUT  POINT,  FRAM  IN  COMMON,  ALSO  AREA,  SIDE  WITH  INTEGRATING  FACTORS 
C 

AREA  *  AREA*(FlOAT(LNOD2)-5.6)/2.4 
SIDE  «  SlDF*FLOAT<LNOD2-4>/4.0 
DO  98  I  »  1,3 

POINT(I)  *  XYZDD(I»11  ♦  ELXYZT(9,I> 

DO  98  J  =  1,3 

98  FRAM(I,J)  a  FRAME(I,J> 

RETURN 

C 

C****  WRITE  DIAGNOSTIC  ERROR  MESSAGE . 

C 

99  WR I TE (6 , 699 )  NEkROR 

699  FORMAT ( /6H  ERROR, 15, 18H  IN  SEGMENT  HALCOF) 

STOP 

END 

SUBROUTINE  VECTOR<U,  V,  U) 

C 

C****  TO  COMPUTE  VECTOR  PRODUCT  U*V  INTO  AREA  W. 

C 

DIMENSION  U<3) ,  V ( 3) .  W<3> 

X  *  3 

DO  2  I  ■  1 ,3 
K3*6-I-K 

U<K3)«U(K)*V(I)-U<I)*V<K) 

2  K  “  I 
RETURN 
END 

SUBROUTINE  SCAlAR(U,  V,  PROD) 

C 

C****  TO  COMPUTE  SCALAR  PRODUCT  OF  VECTORS  U  AND  V. 

C 

DIMENSION  U ( 3 ) ,  V(3) 

PROD  •  0.0 
DO  2  I  =  1 .3 

2  PROD  *  PROD  ♦  U(I)*V<I) 

RETURN 

END 

SUBROUTINE  SFR(XLOCAL,  WCORN,WLOOF,  NSTaGE) 

C 

C*****  SHAPE  FUNT10N  SUBROUTINE  TO  SERVE  HAOOF . 

C 

DIMENSION  MD ( 4 ) ,  TERMV < 46 ) ,  WCORN(10,3),  WLOOF<10,3),  XLOCAL(2> 
COMMON /FROn/ELST (528) ,ELR<32,2) ,LN0DS(8) ,LNODZ 
COMMON/SHEI L/AREA.  ELXYZT<9,4>,  FRAM<3,3),  POINT(3), 

1  SIDE,  TH I K  »  WSHE  L  < 1 3 , 45) .  XITA<2),  XYZ<241,3) 

COMMON/ COE F/COEF (247) 

EQUIVALENCE  (WSHFL<1,1),  TERMV(1)> 

DATA  MD/8,  43,  90.  171/ 

C 

C****  INITIALIZE  AND  PREPARE  TO  CALCULATE  TERNV  =  POLYNOMIAL  TERMS. 

C 

XI  *  X LOCAL ( 1 ) 

ETA  «  XL0CAL<2> 

IF(ABS(XI) -GT.1 .0.  OR  . ABS ( ETA) . GT . 1 . 0 .  OR 

1  .(LNODZ.EQ.6  AND  .  (X I  .  L'  . 0 . 0 . OR . ETA* (1 . 0-XI -ETA) . LT . 0 . 0) ) > 

2  GO  TO  99 
IA  «  2 


«■»<-><•»  o  o  o  non  non 
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****  CREATE  POLYNOMIAL  TERMS  AND  XI*  ETA  DERIVATIVES. 

TERMVd)  *  0.0 

TERMVI2)  «  1.0 

NZ  *  <LN0DZ+NSTAGE-3>/2 

DO  6  N  *  1 .  NZ 

IAN  *  IA  ♦  N 

N2  *  N  ♦  15 

N3  *  N  ♦  30 

DO  4  J  «  I A , I  AN 

TERMV(J*N)  a  TERMV(J)*XI 

TERMV (J*N2)  a  TERMV<J)*FLOAT(lAN-J> 

4  TERMV  ( J*N3)  a  TERMV(.I-1)*FL0AT<J-IA) 

IA  =  IAN 

6  TERMV < I A  +  N I  a  TERMV ( 1 A-1 ) *ET A 

****  CREATE  SPECIAL  COMBINATIONS  FOR  LOOF  NODES,  ETC. 

DO  8  I  *  8.  38.  15 
I  F  <  LNODZ . EO . 6)  TERMV ( I >  a 

1  2.0*(TERMW<I)-TERMV<I+3) )  ♦  3 . 0* < TERMV ( I +1 ) -TERMV ( 1*2) > 
I  F  <  LNODZ . EQ . 8)  TERMV(I)  a  TERNV(I*2> 

I  F  (  LNODZ .  EQ  .  8)  TERMVH+2)  *  TERMVCI+6) 

8  CONTINUE 

****  USE  TERMV  TO  FIND  WCORN  AND  WLOOF  AND  XJ,  ETA  DERIVATIVES. 

NFOISZ  a  (NSTAGE  +  D/2 
DO  18  NFOIS  a  i ,  NFOTSZ 
NZ  *  <3*LN0DZ)/2  +  NFOIS  -  4 
IF(NZ.NE.IO)  GO  TO  12 
NZ  *  9 

DO  10  I  =  10,40.15 

10  TERMVU)  ■  TERMVU+3)  -  TERMV (1*5) 

1 2  K  s  0 

DO  18  I  *  1,3 
DO  16  N  =  1 ,  NZ 
GASH  a  0.0 

LNODZ6«LNODZ+NFOlS-6 
MDEL*MD(LNODZ6>+N*NZ-15*I 
MA  *  16*1-14 
M2  *  15*1 +N2-1 4 
DO  14  M  «  MA,  MZ 

mmdel*m*mdel 

14  gash»gash*termv<m)*coef<mmdel) 

NK*N*K 

IFCNF0IS.EQ.1)VC0RN<NK, I)»GASH 
NK«N*K 

IF<NFOIS.E0.2)WLOOF(NK, I)=GaSH 
16  CONTINUE 
18  X  *  1 
RETURN 

****  ERROR  DIAGNOSTICS,  IF  POINT  LIES  OUTSIDE  ELEMENT. 


A 


99  WRITE (6 , 61 0)  XI,  ETA 

610  FORMATI/30W  ERROR  11  IN  SEGMENT  SFR,  XI  * , F 1 5 .9 , 3X , 5HETA  »,F15.9> 
STOP 
END 


FINISH 
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LIST  OF  PRINCIPAL  SYMBOLS 

flexural  rigidity  =  Eh^/ 12(1  -  v2) 

Young's  modulus 
h2/ 1 2R2 

shell  thickness 
radius 

components  of  displacement  acting  in  the  x,  y,  z  direction 
components  of  displacement  acting  in  the  6,  $  direction 
component  of  displacement  normal  to  the  shell  surface 
components  of  strain 
components  of  curvature  change 
Poisson's  ratio 

components  of  rotation  in  the  x,  y,  z  directions 

components  of  strain  acting  in  the  0,  $  directions 

spherical  polar  coordinate  measured  in  the  xy  plane, 

spherical  polar  coordinate  measured  from  the  z  axis  to  the 
xy  plane 
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LIST 

PROGRAM(DATA) 

INPUT  1 =C  RO 
0UTPUT2sLP0 
OUTPUT3*CPO 
END 

MASTER  GENnATA 

c****  program  to  generate  x,y,z  coordinates  on  thf.  surface  of  a  sphere 

C»***  GIVEN  THE  CORRESPONDING  THETA  AND  PM  I  COORDINATES. 

C****  RrRADlUS 

C****  NNOoE=TOTaI  NUMbE«  Of  nodes 
REAL  X(241) , Y(241 > , Z  <  241 ) 

R=10.0 

PI=3.14159?6S36 
PI F  =  P1 /I  80 . 0 
READd  #100)NNODE 
DO  50  I  NOD E  =  1 , N MODE 
READd  »101>NODE  ,  THE,  PHI 
I  F  (NODE.NE  .  I  NOoi-:)  WRITE  <2, 200  I  NO[>E,  INODE 
C****  CONVERT  THFTA  AND  PH!  TO  RAMANS 
THE=THE*PIF 
PMl  =  PHI*Plt: 

X(NODEIsR*COS<THE)*SINCPHT) 

Y(NODE)=R*^IN(THE)*S!N(PHl) 

Z(NODE)=R*COS(P;|I) 

50  CONTINUE 

c****  READ  IN  AND  CALCULATE  MIDSIZE  NODES  NOT  ALREADY  CALCULATED 
C****  NODE  LIES  HALFWAY  AuONG  THE  GEODESIC  LINE  BETWEEN  NOA  AND  NOB 
60  READ(1  ,102,END  =  70)NODE,NOA,.l0R 
XFACs<X(NOa)*X(:jOB>  ;  **2 
YFAC  =  <Y(NOA)*Y(-:Od)  )**2 
ZFAC*(Z<N0a>*Z(M0B))**2 
X(NODE)sR*s<)RT  ( X  ?  AC /  (XFaC+YcAC*ZFAC)  I 
Y(N0DE)sR*S(5RT(YFAC/<XFACi-V(-  AC-*'ZFAC)  ) 

Z(NOOE)=R*rORT(ZFAC/f.XFAC  +  YFAC  +  ZFAC>> 

GOTO  60 
70  CONTINUE 

00  80  IN0DE  =  1  .  NriODE 

UP  ITE  (3,300)  TNO  ,F.  X  '  »NODF>  ,  vdNODE>  ,  Z  (INODE) 

8't  CONTINUE 
STOP 

10D  FORMAT (10) 

101  FORMAT ( 10, 3F0.0) 

102  F0RMAT(3I0) 

?00  FORMAT  (1H0  ,  'NOD? I^  .»  a  NOD.', I*,) 

TOO  FORMAT (I4,3f15. 1  0  > 

END 

FINISH 


Fig  1  cont'd 


LIST 

PROGRAM<DATA> 

INPUT1*CR0 
OUTPUT3«CPn 
OUTPUT  4*CP1 

end 

MASTER  GENoaTCH 

c****  PROGRAM  TO  GENERATE  X,V,Z  C'iOKDI  NATES  FOR  A  3*3  PATCH  OF 
C****  RECTANGULAR  ELEMENTS  ON  A  SPHERICAL  SURFACE. 

O***  THE  INPUT  nATA  1$ 

C***+  I  TYPE  ICASE  T»EINC  pH  1 1  *jr  THEZFpO  pHTZES-iO 

C****  ITYPE«1  FOR  SEMUOOF, 

C****  ITYPE*2  FOf  NASTRAN. 

C****  ICASEaO  FOR  PRESSURE  LOADING 

C  ***+  ICASE«1  FOP  KThupHJ 

C  **»*  I CASE*2  FOR  aThETHE.KPHIPHI 

O***  1CASE  =  3  FOR  ETH.'PHI 
C****  I C ASE*4  FUR  ETHcTHE, f PHIPHI 

C****  THEINC»DEGRFES  SUBTsNOFD  BY  F IF ME NT  IN  THETA  DIRECTION, 

C**«*  PHI  ^"DEGREES  j'.JBTc'lOEO  BY  ELEMENT  IN  PHI  DIRECTION. 

c****  thezero=thbta  coordinate  at  centre  of  patch. 

o***  PHlzEROsPH!  COORDINATE  at  C-NTRe  OF  PATCH, 
f***.  RsRADIUS  OF  PATCH. 

C***  +  ARRAY  N  CONTAIN.;  SEI1ILOOF  N OnE  NUMBERS  ^POUND  BOUNDARY 

C  * ** *  ARRAY  M  CONTAINS  S2MUOOF  M  0$ME  NODE  NUMBERS  AROUND  BOUNDARY 

c****  ARRAY  NaST  CONTAINS  NASTRAN  NODE  NUMBERS  AROUND  BOUNDARY 
C****  ARRAY  NASTSfnE  CONTAINS  THE  S I  Dc  NUMBERS  FO?  T  H  F  NODES  IN  NAST 
C****  A  CORNER  OF  THE  NaSTPaN  PaT  ;h  HaS  SIDE  =  S 
REAL  X  ( 40) , Y ( 49 ) .2(40) *  A THE < 40) » A  PH  I  < 401 
INTEGER  N  <  1  2)  *  M  ( 1  2)  ,  NAST  (12>  .  NASTSI0FO2) 

DATA  N/ 1,3, 5, 7, 12, 13, 23, 29, 34, 36, 38, 40/ 

DATA  M/2, 4, 6, 11, 22, 3 t, 3°, 37, 35, 30,1°, 8/ 

DATA  NAST /I ,2,3,4, 5, A, 9,12,13, 14, 15, 16/ 

DATA  NASTStdE/5,1 ,1 ,5, 4,2, 4, 2, 5,3, 3, 5/ 
psIO.O 

P!»3. 1415926536 
PIF*PI/180.0 

RFAD(1 #101 > I  TYPE, I  CASE, THE!  iC , PHI  INC , THE ZERO , PH  t  ZERO 
C****  CONVERT  TO  RaDIaNS 
THEINC«THETNC*PIF 
PHI INC*PHI TNC*P’ F 

thezero«thezfro*pif 

PHIZEROsPHT  ZF«0*PI F 

C  *  **♦ 

GOT  0(0,40) , ITYpf 


C  *  **  * 


TM  Str  954/2 


Fig  1  cont'd 


c****  CALCULATE  COORD  T  NATHS  FOR  S.MIL.IOF  PATCH 
PH  I  =  PH  I ZE RO  +  1 . 5*PHI INC 
NODE*1 
DO  30  1*1, A 

the*thezero-i.5*theinc 

00  10  J=1 ,7 

X(NODE)*R*rOS(THE)*S»N(PH!) 

Y(NODE)=R** TN(THP)*STN(PHI> 

Z(NODE)*R*'-OS(PH1) 

ATHE(NODE)=THt/p!F 
APH! (NODE)sPHI / p I F 
NODE=NODE*1 
THE=THE+THFlNC/2.0 
10  CONTINUE 

I F ( I . EQ . 4 ) GOTO  <0 
THE  =  THEZER0-1 . 5  *THE i NC 
PpI=PHI-0. 5*PHI INC 
DO  20  J  3 1  ,  4 

X(NODE)*R*m$<T  I  =  >*;)»  n<  PHD 
Y(NODEJ*R*^IN(T;jF)*STN(PHn 
7(N0DE)*R*COS<P.U  > 

ATHE(NODE)=THE/jIF 
APHI (NODE)*P«I/p|F 
N0DE=N0nE*1 

the=the+th=inc 
20  CONTINUE 

PHl=PHI-0.5*pHI INC 

30  continue 

WRITE  (3, 300)  (  NO.)*,  X;  NODE)  .Y^NODE)  ,Z(NODE)  ,NODE  =  1  ,40) 
URJTE(4,400) ICASF 

C****  PRINT  NODE  NUMij -R  ,  THETA ,  PH  I  FOR  EACH  BOUNDARY  NODE. 

DP  32  1=1,12 
NODE*N(I ) 

32  WRITE <4, 400) NODE, AT  |F(NODE>  . A  PH  I (NODE) 

UR ITE ( 4, 401 ) 

C****  PRINT  NODE, THETA, PHi .SIDE  NUMBER  FOR  EACH  BOUNDARY  MIDSIDF  NODE 
<  =  1 

DO  34  1=1,4 
DO  34  J*1  ,  3 
MI D  =  M<  K) 

WPITE(4,4U0)MI D. ATM J  ail D) , A  >H I (  1  \ D ) , I 
3u  Y=K>1 
STOP 


1 


***  CALCULATE  COORDINATES  FOR  N„STRan  patch 
4  PHl  =  PNIZERO+1 . 5*CH I  INC 

node31 

pO  60  1*1, 4 

THE  =  THEZER0-1 . 5*the  x  m c 
pO  50  J  =  1  ,  4 

X (NODE) =R*rOR (THE) *  3TN(PHT  > 

Y(N00E)aR*SlN(T-Ifc')*5TN(PHl> 

Z(NODE)*p*COR(Pf)I ) 

athe</vod€)sthe/;>i  f 

ApHl  (NOOE)'spHI/.'Ic 
N  0  D  E  *  N  0  D  6  ♦  1 
the=the*thfinc 

SO  CONTINUE 

°Hl=PHl-PHT INC 
60  CONTINUE 

*.*  PRINT  NODE  X  Y  Z  AT  FaCH  NO  )E 

WRITE  <3, 30?)  ( INODE  ,  A<  INODE)  V  <  I  *|ODF  >  ,  Z  <  T  NODE  >  ,  INPDEsI  ,  1  6  > 

**.  PRINT  ICASc 

WRITE(4,40?) ICAO* 

».»  PRINT  NODE, THETA. PHi  AND  S I  > F  NUMBER  FOR  EACH  BOUNDARY  NODE. 
DO  70  1=1,1? 

NODE*N4ST<  T ) 

\ ? IDE  =  NASTs  rpE (  I ) 

WPlTE(4,400)N0D,:,AT:Ie(NODF).APHT(N0DE),!SlDE 
7  CONTINUE 
STOP 

1 7  1  F 0 RMAT (2  I  0 , 4P 0  .  )) 

7  ».-■  FORMAT  ( 14,  SF15  .  ’  0) 

7)7  FORMAT  <  'OR!  D'  ,4.<.  13,  AX,  3F8.  .) 

FORMAT  (|4,2F10.,?  ,  I  4  > 

4 oi  format ( •  o  o  o  n.V) 

407  FORMAT  ( !  1  ) 

END 

FINISH 
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C  *  *  *  * 
C**** 
C  *  *  *  * 

c  *  *  *  * 
c  *•*  *  * 
c  *  **  * 
c  *  *  *  * 
c  *  **  * 
o*** 
c  *  *  *  * 
c  *  *  *  * 
c  *  *  *  + 
c**** 
c  *  **  * 


c  **** 
c  *  *  *  * 


c  *  *  *  * 


c  *  *  *  * 


LIST 

PROGRAM(DATa) 

I NPUT 1 *C  RO 
OUTPUT3=CP0 
0HTPUT4=CP1 
PND 

master  gentripat 

THIS  PROGRAM  GENERATES  DATA  FOP  A  PATCH  OF  18  TRI ANGULAR 

sfmiloof  eiements  Ol  THE  S'J-PACF  Of  a  sphere. 

THE  INPUT  <vatA  ;3 

ITVPE  ICASf  THE  t  NO  PHIIUC  T!E?FRO  PMIZERO 
WHERE  ITYpp  IS  DUMMY 

ICASP  IS  THE  BOUNDARY  DISPLACEMENT  TYPE  0.1, 2, 3  OR  4 
theinc  is  the  increme  it  per  element  in  the  theta  DIRECTION 
PHIINC  IS  THE  TNCRFME  ,T  PgR  ELEMENT  IN  THK  PHI  DIRECTION 
THEZFRO  ID  THt  THETA  COORDINATE  AT  THE  CENTRE  OF  THE  PATCH 
PHIZFRO  IS  T  H  t  PHI  COORDINATE  AT  THE  CENTRE  OF  THE  PATCH 
NOOE  X  Y  Z  ARE  WRITTEN  TO  CP  F I LE 3  FOR  FACH  NODE  AND  MIDSIDE  NODE 
NODE  THE  PHI  FOr  EACH  NODE  ON  THE  BOUNDARY 
MID  THE  PHI  T  S I ; .  F  FOR  EACH  HtjSIDE  NODE  ON  THE  BOUNDARY 
ARE  WRITTEN  TO  CP  FILE4  FOR  INPUT  TO  THF.  SEMILOOF  PROGRAM 
RFAL  X<49)  ,  y(49>  ,  ZU9)  ,  ATHE  (49)  ,  APHI  (49) 

INTEGER  N(12) #M(12) ,NAST(12 ' 

N  CONTAINS  THE  VpRTeX  NOnES  AROUND  THE  BOUNDARY  OF  THE  SEMILOOF  PATCH 
M  CONTAINS  THE  UIDSinE  NODE'.  AROUND  THE  BOUNDARY  OF  THE  SEMILOOF  PATCH 
DATA  N/ 1,3. 5. 7, 15. 21. 29, 3 5,. ,3,45,47, 49/ 

DATA  M/2,4  ,<$.  14  .  ?6.  <♦?  .  4.1 . 46  •  44, 36, 72  ,#/ 

R=10.0 

PI*3. 1415926536 
t>lF  =  PI/180.0 

READ (1# 101) t TYPE, I  CASE, THF I  JC, PH  I  INC, THE  ZERO, PHI  ZERO 

CONVERT  TO  RADIaNS 

THElNC»THE!NC*PlF 

PHIINC=PHI INC*PIF 

THEZERO=THcTERO+p I F 

PHlZER0sPHtZF90*°IF 

calculate  coordinates  for  s.miloof  patch 


Fig  1  cont'd 


10 


3  0 

C  *  *  *  * 


3 7 


3^ 

10  1 

3  n  -i 

4  )( 

M1 


PHl  =  PHIZER0  +  1 . S*?H1  INC 
N  0  D  E  =  1 
00  30  1=1.? 

THE=THE7ER"-1 . 5*THE  INC 
00  10  J  =1  , 7 

X(N0DE)=R*C/>S  (T  !F)*  i  IM  /  pm  ) 

V(  NODE)  =R*S  IN  (T  :;c)  *S  ’  M  (PH  P 
7.  <  NODE  )  =R*r  OR  (  p  >  I ) 

ATHE(NODE)=tHE/,'IF 
A PHI (NOoE)-PHl / P  !  F 
NOOF=NODE+1 
THEsTHE+THFINC/’.O 

continuf 

THE=THE7ERO-1 . 5*THPI^C 
PHIsPHI-0,R*dH| ?mc 
CONTINUE 
PRINT  RESUlTS 

UR  ITE  <  3 .3001  (ijOOF  .  X  (MOPE)  ,  Y  •  NODE)  ,  z  (NOOf)  ,  NOOE  = 

WRITE(4,A00)  T  C  a  -)  f 
DO  3 £  1  =  1,1? 

NODE=N<  I  ) 

WP  ITE  U,401)NOfK  ,  AT,(r  CJO0E1  .  A->HT  (NOPE) 
WRIT£<4,401  ) 

X  =  1 

on  34  l  =  1  ,  a 
00  34  J  =  1  ,  7 
MID=M<<) 

W  R  I  T  E  <  4 , 4  0  o  >  M  I  r>  .  A  T  M  E  1 1 1 .0  )  ,  A  H  1  (  1 1  i)  >  ,  I 

K  =  K+1 

STOP 

FORMAT ( 2 1 9 , 4  f  0  .  )  1 
FORMAT  < 14, 3Fl  5  .  0  ) 

FORMAT ( 14, ?p1U.  \  14; 

F0RMAT(*0  O.u 
END 

FINISH 


,49) 
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LIST 

prOGRAM(DATA) 

1NPUT1=CRO 
OUTPUT  3  =  CP0 
0UTPUT4=CP1 
END 

MASTER  GENipaT 

c**.*  THIS  PROGRAM  CAt.FUI.ATES  THE  XV  Z  COORDINATES  FOR  A  3  BY  3  PATCH 
C****  op  SEMI  1.00  c  OR  MaSThaN  OU  A  Dk  '  LA?  E  R  A  L  ELEMENTS 

C»***  WITH  A  r  E  G  (i  L  A  R  ;OijN;)ArY  R'.it  ulTH  IrrEGijI.AR  ELEMENTS  IN  THE  INtErIOr. 
C****  THE  INPUT  DATA  *S 

C*»**  ITYPE  iCASt  T  H  F  j  N  C  ,’HIINC  T  i  E  Z  E  R  0  PNtZERO 

C  *  *  *  *  THEI  PHII 

C  +  ***  (1=14.16,25.27)  '-OR  SEMii.oOF  OR  <1=6,7,10.11)  FOR  NASTRAN 

C*o*.  WHERE  ITYPf  IS  1  FOr  SENILOIF  Or  2  EUR  HaSTRAN  DATA 
C  » *  *  *  I  CAS?  IS  the  DOtlNDARY  DISPLACEMENT  TYPE  0.1,  2, 3  OR  4 

C****  THEInC  IS  THE  ANGLF  S'lRTENDED  PtR  ELEMENT  IN  THE  THFTA  DIRECTION 

C****  P  H  I  I N  C  IS  THE  aNv'LF  S''rTE'JDED  PER  ELEMENT  IN  THE  PHI  DIRECTION 

C  *  *  *  *  THEZprO  IS  TS?  THETA  NDORDTNATE  AT  THE  CENTRE  0 F  ’HE  PATCH 

C****  PHIZPRO  r  ;S  T  H  t  PHI  COORDINATE  AT  ^HE  CENTRE  OF  THE  PATCH 

C****  THEI  and  PHII  ARE  THE  THFTA  AND  PHI  COORDINATES  OF  The  VERtEX 

C**»*  NODES  AROUND  The  CENTRAL  ELEMENT. 

REAL  X  (40>  .  YI40)  ,Z  I.f0)  ,  ATH5- 40)  ,  APHJ  UO) 

INTEGER  N(12)  ,M(l2)  ,  II I  D  (  3f>)  NAST(12).fJASTSl:'E(1?) 

DATA  N/1 . 3, 5 ,7, i2» 1 3 . 23, 2* , 7  4, 36, 38 , 40/ 

DATA  M / 2, 4, 6. 1 1 . 22,. 5  3, 30. 3  7. 3  5,30, 19,8/ 

DATA  MI0/9,3,U.l0,S. 16,13,  .2,14,15,14.  1  4 . 1  ’ ,  1  6 . 1  3 , 7  G  .  14,25. 

1  21, 16, 27,  ?4,23,?e,  26, 25, 27. 28, 27, 20. 31  .25,  <6. 3? , ??. 3A/ 

OATA  NAST/1  ,2,3.4,5,8,0,12,  .3,14,15,16/ 

DATA  NASTSTDE/5,1,1,5,4,2,4-2,5,3,3,6/ 

3=10.0 

Pl=3. 14159?6636 
PIF«PI/180.0 

READ(1,101)TTYpp,IC. 1 SE, THEINC, PHII  NC,tHFZER*i, PHIZ  FRO 
c  *  +  *  *  calculate  coordinates  fop  s  mu  oof  patch 
c  *  ★  *  ★ 


GO  TO  <0 , 75  )  ,  f  t  Y 0  ,E 


Fig  1  cont'd 


phI=pH!2ERo+i  .*5»pHT  ttic 
?J  n  D  £  =  1 

00  30  1=1.  L 

Th£=THE?ER0-1  . 5  *THE  ;  fJC 
no  10  J=1 ,7 
AT«E (NOPE) =  THE 
4  PH  I ( NOPE ) - PK I 
N0DE=N0DE+1 
THE=THE+THr INC/?  .  1 
in  continue 

I f ( I . EQ.Ajnnro  jo 
th£=THEZER0-1 . 5  *THE I NC 
PHl=PHI-0.5*PH1 I nc 
DC  20  J  =  1  .  4 
4THE (NODE) =THE 
,4PHl(H0n£)=PHl 
‘IOOE  =  NOnE+l 
THErTHE+THFlNC 
?rt  rPNTINUF 

OHIrPHI-0.i*DHI  INC 
30  COMTINtJF 

C****  READ  IN  THE  TH  =  ta  ft, In  p;|!  rod  PM  M  aT  £3  Op  The  N  0  o  E  $  OF  THE  CFNTP4L  ELEMENT 
REA0(1.i02)ATHE(i4),APm<14- 

RFA0( 1 , 1 02 ) ATHE( 16) , A ?HI ( 1  4  ’ 

REA0(1.1O2UTHe<25),.»PHI(2S. 

3  E  A  D  (1  , 1  0  2 )  a  T  H  P  (  2  7  )  ,  A  P  H  I  (  ’7  > 

D ft  40  1=1  ,-*4,  r> 

4THE(MIp(I))  =  (athe(;Mo(I  +  11'^aThE"i:d(I+2>))/2,0 
CO  4PHI  (MI  D(  I  )  )  -  {  4  PH  I  (;l»  D{  1  + 1  )  ♦APR  I  (M  I  0(T  *2  >  )  )  /2  .  ft 
no  sO  1=1,40 

y  M  )  =R*COS  <  aTHf  < I  )  *j>  r  F  )  *3  {  N  ■'  a-'H  i  '  I)  *»  :  F  ) 

V  (I)=R  +  SINIATHE  <  I  >  * ,»  If)  *SIN  A<-'Hi<P*PlF> 

3'  2 ( I)=R*COS(APHf (I)*?Tf) 

C  *  *  *  •  n  E  I  NT  R  E  S  U  l.  T  <5 

WP  I  TE  <  5 , 300  )  <Uu  jF  .  X  (MODE)  ,  V  NOnn  )  ,  i  r  rgon1: )  ,  Nnftf  =  1  ,  4  0  ) 

UPIT£(4,409) f  C  A  5  E 

DO  60  1=1,12 

N  ft  0  £  =  N  (  I ) 

6:  WP I T£ (4, 400) NODE , AT  IF (NODE)  .APHi  (NO  OF) 

HP  I TEi4.40l ) 
x- 1 

DO  70  I  a  1,4 
-»n  70  J  =  1  ,  3 
H  I  =  M  C  K  ) 

WRITE  (4,^00)  MI  .  *TrtE.;MI)  ,  APH  (MI)  ,  I 
T  '  V  =  K  ♦  1 
"TOP 


Fig  1  cont'd 


o***  CALCULATE  ''OnRD;\ATc$  FOR  N.,STRAN  PAtCH 
75  PH]=PHUER0  +  1  .5*PH!  INC 
N0DE=1 

oo  as  1=1,4 

THE=THE?ER0-1 ,5*THE,NC 
pO  80  J  = 1 , 4 

athe (NOpe) =the 

APHI ( MOPE ) =PH I 
N0DE  =  N0pEC 
THE  =  THE  +  Trt<:  INC 
80  CONTINUE 

PHI=P«I-PH7 INC 
8S  CONTINUE 

R  F  A  D  ( 1  ,102)  A  THE  i  6  ^  » .  '« p  H  I  (6) 

R E A 0 ( 1 ,102) 4  T  M  f  7) ,  n  H  I ( 7 ) 
RFA0<1,102)ATME('<C.  >,ApHlf10 
READ<1 1 102) A* hf  /  1 1 )  ,  ApHI  (1 1  . 

DO  90  1=1,1? 

J=NAST(I) 

9  0  UR  ITE<4, 400  I. >, AT  Ht(J>,APMl(.->,NASTSIPE(I) 
c****  CALCULATE  aNP  PRINT  the  X,Y.Z  COORpINaTeS 
DO  95  1=1  ,  1  A 

X< I ) *R*COS  < ATHE ( I >*pf  F>  *S I N  .  A°HI  '  I  >  *P  I  F  > 
Y(I)=R*SIN('ATHFfI>'»PlF>*SIN  APMi<I)*Pin 
7.(1  )  =8*C0S  (  APHI  (I  >  *PT  F) 

9  5  WRITE (3, 302) I , X(  I ) , Y(  I ) , ?( I* 

STOP 

101  FORMAT (2I0,4f0.0> 

10?  FORMAT  (?F0  . ;)) 

30P  FORMAT< I4.3F15 . 1 0> 

30?  FORMAT ( • GR T p  *  »  4X , I  8 , 8X . 3F8 . 0 ) 

400  FORMAT  <I4,2f1Q..?, 14) 

401  FORMAT  CO  0.0  J.OM 
END 

FINISH 


Fig  }  cont'd 


LIST 

PROGRAM(NAST) 

INPUT1=C»<» 

OL'TPKT3  =  Cpf> 

O'JT  PUT4*C  P 1 
OUT  PUT5*CP? 

EN  D 

MASTER  0  E  K V  * S  T  a  ; 

c.***  this  program  tip  :rkat  =  s  the  ioundary  displacements  for  a 

c****  NASTRAN  PATCH  .11/  ,iY ,  UZi’IX  ,'J'  ,’tz  ARE  PRESCRIBED  AT  EACH  NODE 

c****  on  the  boundary  n  the  form  or  sec  hr  m»c  cards. 

C  *  *  *  *  ICASE  is  RF  AD  I  J  A  "I  i  THE  APPROPRIATE  SUBROUTINE  13  CALLED. 

C****  JCASE  =  0  P  R  F  S  S  U  3  F.  L'i.^IMO 

€♦♦**  !  C  A  S  E  =  1  I  N  c  Y  T  F  h  i  I  0  H  ■  OJNMifi.  y.APPA  PH  f  T  H  E  r  A  CONSTANT 

c « *  *  *  i  c  a  s  e  =  ^  infxten  ;  r-)N,v  rcNoi  !<?#*, \ppa  phi  &hi, -kappa  the  the  constant 
C  >  *  *  +  I  C  A  S  E  =  3  MEM", RAN-:  iic  lAVJU'.IR  .  ;£>?  I  I.GN  PHI  TMETa  constant 

C»***  ICASE=4  MENRPAN-  jRHAVIiMIP,  :pSl(.ON  PHI  oh  I , -EPS ' LON  THE  THE  CONSTANT 
f>*  *  *  P  =  R  APIUS 
C****  h  =  TH  ICKNESS 
C  «  »  *  -  “^YOUNGS  MODULit  > 

c  *  *  +  *  phu=poi ssoks  Rarrv 

c.***  NOD  =  NUMAE  H  OF  NODES  On  •;  O  U  N .  i  A  -H  Y 
c**»*  irrtotal  number  of  (Odes  in  patch 

C  *  *  *  *  PIF  =  PI/180  =  MULTIPLICATION  c  A  C  t  n  r  to  CONVERT  FROM  DEGRESS  to  RADIANS 
C  *  *  *  * 

R  =  1  O  .  0 

H=0.04 
?=1 0 . OEft 
DNU=0.3 
‘JODal  i 
IP  =  16 

PI  f  =  3 . 141  5u?65Ti/ i  )  1 
rFAD  < 1 #100)  ICAs? 

I  F  <  1C  ASE  .  ED  .  o)  r  LL  .'RcSSlIKc  R,NOD»pIF) 

I  F  (  ICASE.  EO.  1  >Csl  n  THE  PH  l  C  •  ,  H  .  c  .  PM"  r  NOn  ,  I  P  ,  P  I  F  ) 

IF  ( ICASE.  t0.2>C  ALl  r.HFTHE  (  » ,  NOd  ,  I  o  ,  e  I  F  ) 

T  F  (  ICASE  .  60. 3>  C.'.LL  JTHEPHI<--.NOd,Ip,pIF) 

I  F  (ICASE.  EG  ,4)C;,LL  ;thcthE<  : .  H 0 f.  ,  I  P  ,  o  I  c  ) 

STOP 

1),,  C0RMATO1) 

END 

subroutine  pres  ;ur=  i-'.nop, p  o 


Fic  1  cont'd 


C***  +  ICASE  =  0 

c»***  THIS  SUBROUTINE  GENERATES  T  if  BOUNDARY  CONDITIONS  for  a  3  by  j  patch 
C****  OF  NASTRAN  0IIA04  ELEMENTS  U  'HER  PRESSURE  LOADING 
C****  THE  input  DATA  can  JF  OfiT  A  t  ;n  FROM  PROGRAM  GENFATCH 
C**»*  THE  OUTPUT  IS  UX,‘JY,IIZ,:IX#U'VJZ  pO*  EACH  NODE  AND  IS  WRITTEN 
C****  TO  CP  FILE'?  AS  |oC  CARDS 
DP  50  11*1. NOD 
READO  .  1U1M  .T4E.PHI  ,  IS  IDE 

C***+  1 S I P E * 5  REFERS  TO  A  NODE  ON  THE  CORNER  OF  THE  PATCH 
THE*THE*P I c 
DHI*PH1*PI  F 

GO  TO<13#24 ,1 5 .24. 13' , I  SI  of 
13  CONTINUE 
C  *  *  *  *  SIDES  1  AND  3 

AX*”COS(PHI ) *  C  0  S ( T  H  e  Y 
AV*-COS(PHn*Sl  JYTHE) 

AZ*SIN<pHO 
rx=-sin(Th?y 
RY*COS (THEY 
RZ*0. 

WRITE (3, 302) T , A”. I , AV,  I 
WRITE(3430?)I,I.sZ 
UR  I TE ( 3 » 300 ) ! .Ha. I »  RY  ,  I 
UPIT£( j.301 M  .  F , RZ 
IF<  ISIDE.EO.  SJGiiTU  24 
r,n  TO  50 
24  continue 
C  *  *  *  ■*  SIDES  2  AN  *  4 
AX*-SIN(THF) 

AY«COS<THF) 


AZ*0  . 

PX*-C0S<PHI  )*CO',(THg) 
RY*-COS(PHT )*SI  J(THC) 
RZsSIN(PHI) 

WRITE (3. 304)  I  .A''.  I  .AX.  I 
WRITE (3. 30S) T , I  ,  AZ 
WRITE  <3,30(*)  t  ,  R  v  ,  J  .  RX  ,  I 
WR I Tg < 3 , 307 ) I , I  .  wZ 
50  CONTINUE 

101  eORMAT(IO,2fO.O.  10) 


300  FORMAT ( • MPC  1 

•  ,  1-3  .  '  4 

1 

,  F  8 . 6  ,  I  8  ,  ■  5 

• , F8 . 6 

18X, ' AMP • , I?) 

301  FORMATCKMp*  ,  I2.11X,  13. 

'6 

*  ,  F  3 . 6  ) 

302  FORMAT(*MPC  1 

•  ,  I  H  •  '  1 

1 

,  F  8 . 6 , 1 8  ,  •  2 

' , F8.6 

18X. 'AMA* , 12) 

30 t  FORMAT(«AM*« ,12f1lx,»s# 

•  3 

'  ,  F  H  .  6 ) 

304  F«RMAT('M PC  1 

•  ,  I  /. .  1  ? 

1 

.  F8.6, i8, ' 1 

' , F8.6 

1 8X . ’ 8MB ' » I 2 ) 

305  F0RMAT<*»MB* , 12,1 IX, t8. 

•3 

'  .  F  8 . 6  ) 

306  FORMAT  < ' MPC  1 

•  ,  I ...  *5 

1 

rFft.6fI8/*4 

' . F8.6 

18X*  '*RC  ,  12) 

30?  FORMAT!  *«MC  ».  12, 1  IX,  TA,  »A  '.F,H.o> 


STOP 

FND 


Fig  1  cont'd 


SUBROUTINE  KTHRPH I ( R . H , E , PNM , MOo , I P , P I  F > 

C***+  ICASE  *  1 

C*«**  THIS  SUBROUTINE  GENERATES  NASTRAN  ^OUNDaRT  DISPLACEMENTS  FOR  THE 
C***  +  INEXTENSIONAL  BENDIIJG  C.\*E  <APPA  CONSTANT 

C****  THETA  PHI 

REAL  U ( 3) 

B  =  R*R*<1.  +  PNII)/<E*H**3) 

DO  10  I a1 *  NOD 
READ<1»110)N.THE,PHi 
TMEsTHE+PI e 
phi=phi*pif 
CPHIsCOS(PHI) 

SPHI=SIN(PHI) 

CTHE=COS(THE) 

STHE»SIN!THE> 

tl(1)=4.0*B*CPHI*CTHc*CTHE*C'HE/SPHI 
U(2)»-4.0*b*CPHI*STHE*STHF*';,  THE/SPHI 

U(3)»-B*COS<2.0  +  THE)*<SPHI*-.PH  i-2.0*CPM!*CPHl)/<SPHl*SPHI> 

DO  5  J*1,3 
if  =  1  00  *  J  +N 

WRITE(3,300>N. J ,U< J) , K 
5  UR!TE<3*301)K 
IP»I?+1 

WPHI*-4.0*B*COS(2.0*THE)/(R»SPHI*SpH!*SPHI) 

WX*-STHE 

WY»CTHE 

UR  I T  E ( 4 , 40S ) N  <  WX  »  IP 
WRITE(4,41S) IP, M,WY, Tp 
UR  I TE! 4, 425) I P. I P.UPHI , IP 
URITE<4,435) TP 
WR I TE ( 5  #  500 )  TP 
I PaI P+1 

UTHE  =  -4.0*B*CPHI  +  SIN  <  2 . 0*TH.')  /  (R*SPM!*S3H  I  *SPH  I  ) 

wx*-cthe+cphi 

wys-cphi*sthe 

W  Z  *  S  P  H  I 

WRIT6(4,400)N.WY, IP 
WRITE(4,410)  TP,;|,UX,  TP 
WRlTE(4/420)  TP . M,WZ, TP 
W0ITE(4,430) TP, TP.WTHE 
WRlTE(b,500)  IP 
10  CONTINUE 
STOP 

110  FORMATUO,?FO.O) 

300  FORMATCSPC*  • #15X, M • ,2'14. F16.9, 'gSPC' , 13) 

301  FORMAT! • *SPC ' , 13) 

4  00  FORMAT! 'MPC*'  ,1C-X,  ’  1  ’  ,  1 1 6 , 1  5X  ,  ’  5  •  ,  Fl  6 .9  .  '  RMPC1  •  ,  I  2 ) 

41 n  FORMAT! *+M PCI ' , 1 2 » 1 X , 1 1 o , 1 5 v . • 4  * . E 1 6 . 9 , 1 6X , ' AM PC  2  * , I  2 ) 

A05  FORMAT  !•  MPC*  '  »10X,  '1',I16,1e:x,,4,»E16.9,'&MPC1,,I2) 

41*  FORMAT!  •♦MPCI  ' , I  2 . 1 X . 1 1  o , 1  V; , • 5 ' , E16 .9 , 1  fix  * ' (MPC2 • , »2> 

4  20  FORMAT!  '*Moc2*  ,  I?,‘I'*X,I1<S.1SX,  *6*  ,  E 1  o  .  9  ,  •  fcMPC  3  •  ,12) 

4  2*  FORMAT! '+MPC  2'  .  1 2 , 1 ?X ,  1 1  6 ,  1  .,x  ,  E 1  A .  9  ,  '  AM  PC  3  *  ,  I  2  > 

430  FORMAT <'*MoC3' ,t?, IX. 1 16, 16;, El  0.9) 

43r,  FORMAT!  •♦MocT' ,  i2) 

50  a  FORMAT! 'SPC* ,12<, ' 1 • , 1 8 , 1 3x  -  •  1  .  .V  ) 

EMO 


Fig  1  cont'd 


SUBPOUT  I  NE  KTHfTH£((<.NOi>,  !  P  >  P  l  F  ) 

C***  *  ICASE  *  2 

C****  THIS  SUBROUTINE  GENERATES  N,5TR*n  ROIJND.hRw  DISPLACEMENTS  FOR  THE 
C****  INEXTENSIOhAI.  RENDllft  CaSE  CAPPA  ,-KAPPA  CONSTANT 

C****  THETA  THETA  PHI  PHI 

REAL  U<3) 

DO  10  1=1, NOO 
R  E  A  D ( 1 ,110>N, THF.HHi 

the*the*pif 

PHl»PHI*PlF 

cphi=cos<phi  ) 

S0HI=SIN(PWI) 

CTHE=COS (THE) 

STHE=SIN(Tme) 

Ul1)s“CTHE**T*(i't>H!*CpMI+i.ii)/<5.0*SPHI> 

U(2)*STHE**X*(Cr'Hl*ijPHI*1.0;/(3.0*SpHl) 

U(3>sCOS<2.0*ThE>  *Cph!*<1 . 0-1 , 0/ ( $PH I *Sph I >  > / 3 . 0 
DO  5  J  =  1*3 
K  =  1 00* J  *N 

WrITE(3,301)n. J ,U(J)  .K 
s  W»ITE<3,301>« 

losip+i 

WPHI**CPHI*COS(2.0*7HE)*(1  .  \  +  !SPHI*SPHl))/(3.0*3*SPHI) 

MXa-STHF 

WYaCTHE 

WRITE<4,A05)N*WX. IP 
URITE(4(41S) tP(M(WV# IP 
WRlTE(4,425)!P,tP,VPHI,lP 
WRITE<4.435) IP 
MPITE(5,50D) tP 
IDSlP+1 

WTHE  =  -SIN<2.0*T.,E)* <  P.O-SPH ; *SPHl-4.0/(SPHI *SPHI ) ) / (6.0*R*SpHI ) 

wx=-cthe*cphi 

UY=-CPHI*STHF 
U2*  S  PH  I 

WRITE<4,400)N.WY,  IP 
WRITE(4,41  0)  IP,.»,yX,  rp 
W9ITE(4,420)  IP..J.MZ,  IP 
MR  I T  E  <  4 , 430 )  T  P ,  i  p  *  ’«ITHE 
WRlTE<5,50D)  IP 
10  CONTINUE 
STOP 

110  FORMAT(TO,?eO.O) 

300  FORMaT('SPC*  ' ,1pX, *1 • ,2 . 16,e16.V, '*SPC’ ,  13) 

301  FORMAT! •*Soc' . 13) 

4  00  FORMAT  (  »MPf  *•  ,1*>x,  •  1  * ,  !16,1.'IX,  •  S'  ,e16.9  .  'JiMPCl  '  ,  12) 

410  FORMAT!  «*MPCi  '  ,  t*.1X  -  I1'i*1  V;*  *4’  , 16X»  'RMPC2'  ,  12) 

405  FORMAT  ( 'MPC*'  #1<>x,  *1  •  ,  I1A,1r.X,  *4'  ,  El  ft. 9.  '&MPC1  '  ,  12) 

41 S  FORMAT ( **MPC1 * , 1 2 , IX . 1 1 6 , 1 5x , • 5 • , El  6 . 9 , i XX , ' *MPC2 ' , I  2 ) 

420  FORMAT < «*MPC 2' ,  |  2 , 1  ?X ,  1 1  A,  1  5X  ,  •  ft  •  ,  j:  1  ft .9  ,  »  *MPC 3  •  ,  1 2 ) 

4  25  FORMAT!  '*MPC  2'  #  (  2 » 1  ?X  » 1 1 ft , 1  .ix  »  E I  A .  0  ,  '  SMPC3  '  ,  I  2 ) 

4  30  FORMAT!  '*MPC3'  ,  I  2 , 1 X  .  1 1  ri,  1  AX  .  El  ft .  9  > 

435  FORMAT! **MPC3' , I?) 

SOD  FORMAT! • SPf • , 12X, ' 1 • , 18,1 3X ■ ' 1 .0* ) 

END 


Fig  1  cont'd 

SUBROUTINE  ETHe  PH  I  !  R  .  NO«> ,  IO.PIF) 

C  *  *  *  *  I  CASE  a  3 

c**»*  THIS  SUBROUTINE  GENc°ATeS  N,STRAN  BOUNDARY  r,]  SPl.ACEYENTS  for  the 
c  *  *  * ♦  MEMBRANE  SOLUTION  EPSILON  constant 

C  *  *  *  *  T M E T A  P H  I 

REAL  U  <  3 ) 

DO  10  lsl.NOO 

P  E  A  D  ( 1 ,110) N, THE  •  P  H  i 

THE=THE*PIF 

RHlaPHUPIe 

Ct>HI=COS(Pwi> 

SPHI=SIN(PHI  ) 

CTHE»COS (TkE> 

STH£=S lN(TMt ) 

SPM I 2=SPH I +  SPH l 

'J(1>»(-(4.rt-p.o«SPMi?  +  SpHI?-*S»HI2)*CTHE*COSf  2.0*THE)*(2.0*SPHI2) 

1  *  (2.0*<!THE*:',  In(2 . 0*THF)  *S?H  I  *Cthe*C’‘S  (2 .0*T‘<E)  )  )  *CPHI 

2  / <6,U*SPHI2*SPH  I ) 

l>(2)s{-(4.0-2.0*SPHi ,?  +  SPMI2*<!pHI2)'*STHE*COS(2.0*THE)-(2.0*SPHl2) 

1  *(2.0*rTHE*3I  N(2 . 0*The)-Sdh  i*ST>1E*C'1S(2 . 0*T,IE)  )  ^  *CPHI 

?  /!6.0*SpMl?.e.pHl) 

IK  3)  =  <  1 .0  +  0  pH  f  *  ‘.pH  I ;  *C'Jn  <2 .  0*THfc )  f  <  2  .  n  *p,  pH  I  ■>) 

DO  5  J*1  »  3 
<r100*J+N 

WRITE(3,300)N* J ,U!J;  .X 
5  WR I T  E ! 3 , 301  )K 
I P* I P+1 

UPHIa-(-(4.0-2.o*SP.il2*SPHI.!*SPi!2)*SPHJ*C0S(2.,i*TH€) 

1  *3.0*(CPHI*CpHl*1  0>*cos-:  2. 0*t-h  £))/(*. 0*?*SPHI2*$PMI2> 

MX*“S  .'HE 
WY=CTHE 

URITE(4,40S)N.UX«  IP 
MRlTE<4#4lS)»f,#‘i.  JY ,  IP 
WRlTE(4,425)IP#?°»')pHI*IP 
URITE(4,43S>  tp 
urITE(5,500) TP 
I p* I  P  +  1 

WTHE*-CPHI  +  (2.0  +  SPHI?>*(1.')-SPHI>*STN(2..')'*TlE)/f3.0*R*SPHI2*SPHI2) 
WX*-CTHE*CPrl  I 

wy*-cphi*sthe 

VI?  ®S  PH  1 

WRlTE!4,400)N,gv. IP 
WPITE!4,411) IP.n.WX, tp 

WPITE(4,420)  IP.-j.W2,  t  p 
WRITE(4,430) TP. iP.WTHE 
WPITE<5,501) td 
1  •'i  CONTINUE 
STOP 

IV'  FORMAT!  !O,2FO.0j 

SO'1  CORMATCSPC*  •  »15X,  *1  •  ,2  lo.Elo.'?,  'ASPC  ,  13) 

TO1  FORMAT!  •*SPC' . 

Ml'  FORMAT!  'Hpo*'  ,  l"x.  ’  .  '  ,  HA,  1  'X  .  '  S'  .  E14.0  .  'RMpCI  '  .  I  2) 

4  V'  Ff)RMAT(  :  2,1.; .  IV, ,  1  S  ,  •  4*  ,E1.,  .0  .  *,^x  .  *  AMP<  2*  ,  I  ?) 

40*  FORMAT!  'MPr*'  ,  V.'X,  '  j  •  ,  I1<S,  1  OX  .  '  .  FVS  ,<J  ,  '&MPC1  '  ,  12) 

41  S  FORMAT!  '  «MpC1  '  .  2, 1„.  Ho,  1  5; ,  '  S'  ,  E  V,  .0, 1DX,  'RMPC2'  ,  T2) 

4  20  FORMAT!  •♦MPC?’  .  :  P  /  I^X,  I1A,  V  X,  '  S'  .E-iD.O,  'AM-->C3'  .  12) 

4?S  FORMAT! ‘*MPC 2' .  '  > .  W. ,  1 1  A  .  1  >X . E 1 6 .9 , ' NMPC 3 '  . I  2 ) 

430  FORMAT!  '*MPC3'  ,  ,  ?  ,  1 .  1 1 .. ,  1 6  ,  Hi  0.9) 

43*  FORMAT!  '*Mdi;3'  ,  >2) 

SO 'J  FORMAT!  *SP'* '  ,  1?;;,  '  1  '  ,  13, 1  3X  .  •  1  .  )'  ) 

end 


Fig  1  cont'd 


SUBROUTINE  KTHETHE  (*  ,  NOd,  I  P  .  PI  F) 

C»***  ICASE  =  4 

o***  THIS  SUBROIITTNf  GENcRaTE  Na'.TKAN  BOUNDARY  DISPLACEMENTS  for  the 
c»***  MEMBRANE  SOLUTION  EPSIL'JN  ,-EPSILON  CONSTANT 

C****  THETA  TRET/,  PHI  PHI 

REAL  <J<3> 

DO  10  1*1, NOo 
READd  ,110)  N, THr, PHI 
THE*THE*PI  F 
PHI*PHI*PIF 
C  PH  I *CO$ (PHI  ) 

SPHI*SIN(P'-I  T  > 

CTHEaCOSITHFI 

STHE*SIN(Tme'' 

U(1  )*<3.0*«;p«iI*‘',PHI*rTHE*CO';  ( ? . 0*Th E ) -2 . 0*COS ( 3 . 0*TH E ) )  / 

1 (3.0*SPH!**3) 

ij(2)*<3.0*SPHl*‘,  phI-»SjHE*CO',(2,i‘)*tHE)-2.0*SIN(3.0*THE)) 

1  / (3.0*SP4 (*»3) 

U(3)*CPHl*rOS(2.0*T||F)/(SDH[*SPHl) 

DO  5  J  *1  ,  3 
<=1U0*J*N 

UR  I TE ( 3 , 300 )N,J,U(J) . K 
s  WRlTE(3,30l )K 
IP*IP*1 

WPHI»+4.0*CPHI*EOS(2.0*THF) /(3.1*R*SPHI**3) 

UX*-STHE 

wy-cthe 

WRITE(4,40S)N,WX, I? 

URITE<4,41 S) j  p , M , UY , T  p 

UrITE(4,425)!P. I P, UP HI , IP 

WRIT€(4,43S)  IP 

WRITE(5,500) TP 

IP*IP+1 

WTHE*0 . 0 

WX*-CTHE*CPHI 

UYa-CPHI*S-HF 

UZ-SPHI 

URITEI4, 400)N,WY,  IP 
WR I TE ( 4 , 41 0 ) 1 P , N , UX , ip 
URITE(4,420) tp, h,UZ, Tp 
URITE(4,430) IP, IP.UTHE 
UR ITE ( 5 , 500 )  IP 
10  CONTINUE 
STOP 

100  FORMAT ( 1 0) 

110  FORMA T<IO,?Ffl.O) 

300  FORMAT  < • S  PC  *  • , 1 5 X , 1 1  * , 2  . 1 6  ,  E 1 6 . 0  ,  •  * S PC •  ,  I  3 ) 

301  FORMAT  (  **SPC  ’  »  I  *>) 

400  FORMAT  <  •  MPC *  •  , m , • 1  * , 116 , 1 Sx , ’ 5  * . FI  6 .9 , • «MPC1 ' , 1 2) 

410  FORMAT  <  '  *MpC  1  '  ,  i2,1;<,  116,  iSy,  *4.  ,£16.9 , 1  ax,  '  RMPC2'  ,  1 2) 

405  FORMaT ( ' MPr  * ' ,  1  , ' 1 ',I16,15X, '  4  '  »  Fl  6 . 9  ,  'RMPCl '  ,  12) 

41 5  FORMAT ( ' *  M  P  C  1 * , l 2  *  1 X • 1 1 6 , 1 5  x , '5' , E 1 6 . 9 , 1  6  X » '  A  M  p  C  2  '  ,12) 

4  20  FORMAT  <'*MPC?'  ,  i?,17X,H6,1*»X,  '6'  ,  FI  6.9  .  '8MPC3' ,12) 

425  FORMAT  < ' *Hpr? ' , 1 2 , 1 ?X , 1 1 6 . 1  AX , El  6 .9 , • 8Mdc3 • , 1 2 ) 

430  FORMAT  I ' *Mpc3 ' ,l2,1X,I16,16x,Fl6.0) 

435  F0RMAT(,*MoC3' , 12) 

500  FORMAT <• SPf • ,12x, *1 • , 18, 1 3X . • 1 . 0' ) 

END 

FINISH 


Fig  I  cont'd 


LIST 

PR0GRAM< DAT/p 
INPUTlsCRO 
OUTpUT3aCPO 
END 

MASTER  GENO')IJAOa 

10  READ<1.100  ,  EMOsOO)  SI,  Ml  #  *!  2 »  M  ’  ,  H  4 
UR  I TE  <  3  ,  30" )  N  ,  N  !  ,  «.? ,  N.3  ,  ',’4 
GOTO  10 
99  STOP 

100  FORMAT (510^ 

300  FORMaT<  •CQUao<*  »,I.3.7X,*1  •  .448,  •  0.*) 

END 

FINISH 


LIST 

PROCiRAM(OATA) 
inputi =CRO 
OUT  PUT3*C  PO 
END 

MASTER  GENCTPI47 
1  0  R 6  AD  (  1  ,  1  00  .  FNOaOO  )  N  ,  n  1  ,  J ?  ,  N 
URITE(3,30P)N,N'i  ,  M  ^  ,  -M  7 
GOTO  10 
90  STOP 

100  FORMAT  <  4 1 0 1 

300  FORMAT  < • CTr I  A3  • , 1 3 . 7X  * ’ 1 ’  . 3  f  8 , *  0.’) 

END 

FINISH 


LIST 

PROGRAM  <  DAT  A ) 

INPUT1 *CR0 
OUTPUTZsCPO 
OUT  PUT3  =  LPO 

CM  Q 

master  gengrto 

INTEGER  I  B C  (  2 5 0 > 

rFAL  X  (250)  ,  Y  (  2  S  C  >  » 2  <  2  5  0  > 

REAO(1  |102)NNO')i: 

DO  10  Isl/NNOUF 

READ  Vi  .100)N,X(*J)  ,  Y  (.  .M )  ,7.<N> 

I BC  < I ) *0 
10  CONTINUE 

20  REAPd » 1 01 .EN0*25) H, !OC(N) 

GOTO  ZO 

25  DO  30  1*1,223 

I f ( 1 8C ( I) . FO . 0 ) MR l T J ( 2 , 200)  ♦ X ( I > , V < I >  ,  Z (l  > 

30  IF  (IBC(I  )  ,NE.  0)  MR  ITS:  <2*  2011  t  .X(I  )  .  Y(I  )  ,?.(  I  )  * IBC(I  > 
STOP 

100  FORMAT  < 10 , 3F0 . 0 ) 

101  FORMAT  <2101 

102  FORMAT(IO) 

201  FORMATdGRTO  ’  .  1.1.  IXr  3F8.  (S.8X#  13) 

200  FORMAT  (  ' GR I  n  •  ,  4X  ,  U,  8X ,  3F8 . 6) 

END 

FINISH 


Fig  1  concl'd 


1 


i. 

V 


i 


I 


► 


C  *  *  it  * 

c  •  *** 
c  *  *  *  * 

c »  *  *  * 


c »  *  *  * 

1  i, 

C  ‘  *  ** 


c  *  *  *  * 


2-'i 

c »  *  ★  * 


40 


10-' 
1  O'! 
102 
107 
70-"' 


LIST 

PROGRAM  (  DC p;j  > 

TNPUT1=CR0 

OUTPUT3sCPO 

END 

master  change  D.iTA 

program  TO  ALTpt<  TM^  COORDINATES  OF  A  SEMI  LOOP  DATA  FILE  SO  THAT 
POSITION  0 F  THE  MlIj'DE  NOD'S  IS  HAI.EUAY  ALONG  THE  STRAIGHT  LINE 
JOINING  THE  TWO  VE9-PX  ?jO.,r 
REAL  NU,X(?O0) ,7<20j) , Z ( 2 0  > ■ 

integer  e l c a > 

READ  IN  STahT  Of  D  4  T  ft  Flic 
REAOU  ,100>NFL 

R  c  A  D  (  1  #101  )  E  » i'i  l.l ,  T  H  l  ,  p  R  tf  S  S  "J  -!  E  *  0  E  N  S  I  T  V 

R  E A  D  < 1 , 1001  NAM 

READ  IN  X  ,  Y  ,  Z  C<»0-tf>(  NATl* 

DO  10  1=1. NAM 

RPAD(1.102)jfx(  { 1  #  Y  <  T  >  #  •»  n  > 

TOR  EACH  El  EMENT  READ  (  J  ELEMENT  CONNECTIONS 
DO  ZO  TNELsi 
REACH,  103ic» 

IF<EL<«1  .f-a.O)k=f> 

REC  AL-.IM.AT  :  RO  S  ,T  i  •)  4  Op  MID.(»F  NODES  APOJNO  TrfC  ElfMENT 

D  o  2  0  J  a  t ,  *  ,  2 

NV1 =EL< J-1 i 

If  (  J  .  NE.O  V'>;?  =  f  .(  !♦  1 1 

I  F  <  j  .  F  Q .  K  )  •.-  v  Z  8  c  <  '<  > 

N"=EL(J) 

X  (  NM  )  »  <  X  <  N  V  1  1  ♦  <  /  N V  ■> )  '  /  p  .  o 

V(NM)«(V(Nvl>*Y^NV2)i/2.0 
Z(NV)a(2(N,'i>  +  7.  NV2)  >/2.0 
WRITE  NOPE  JiMgtH.X.v.Z  FOR  EACH  NOpE 
DO  40  I  =  1  ,  N  M 

URITE(3, 502)1, X«  D,YM),Z(I> 

STOP 

FORMAT ( 10) 

FORMAT (EU.Q.4F0  1) 

FORMAT  <IO,7fO.0j 
FORMAT (810) 

FORMAT  (110, 7f20  1-j) 

ENO 

FINISH 


*  ♦  *  * 
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